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EXPANSION OF DETERMINANTAL EQUATIONS INTO 
POLYNOMIAL FORM* 


BY 


HAROLD WAYLAND 
California Institute of Technology** 


I. INTRODUCTION 


Determinantal equations of the form 


| + +--- +A,| = 0, 


where the coefficients A; are square matrices of order m, arise in quantum mechanics, 
electrical circuit theory, the theory of small vibrations, and in many other branches 
of physics and engineering. In the process of solving such equations, it is frequently 
desirable to expand them into polynomial form. Since most of the techniques for such 
expansion are not discussed in the standard works on numerical computation, it has 
seemed advisable to make a critical comparison of the methods available. The most 
important methods will be described in sufficient detail to aid the non-professional 
computer. 

1. Basic techniques of solution of determinantal equations. Three basic techniques 
for solution of determinantal equations must be considered: (1) Direct solution of the 
equation by numerical methods [1]. (2) Direct solution by the method of matrix 
multiplication [2-6] (directly applicable only to the case | A —I\| =0). (3) Expansion 
into polynomial form [6-14] and solution of the polynomial equation by standard 
methods. 

In spite of the fact that much effort has been spent in developing techniques to 
avoid expansion into polynomial form, that very technique frequently proves to be 
most economical of effort. Let us consider the numerical solution of the simplest case: 


ay * Cie 
a — don 

DO) =|4-—n| =| ™ = 0. (1) 
Gas ane Onn — 


If we assign a value to d and evaluate D(A) using Chid’s expansion [15] of the numeri- 
cal determinant, we shall make the following number of operations: 
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(1/3)(n? + 2n — 3) M-D (multiplications and divisions), 
(n/6)(m + 1)(2n + 1) A-S (additions and subtractions). 


The number of additions and subtractions includes the ” subtractions required to 
evaluate the diagonal terms of (1). If we have D(A) in polynomial form, we can evalu- 
ate it by synthetic division with only » multiplications and m additions and subtrac- 
tions. We can obtain the polynomial expansion of (1) by Danielewsky’s method 
(§2c) with m*—2n+1 multiplications and divisions and n(m—1)* additions and sub- 
tractions; hence to obtain k different values of D(A) by first obtaining the polynomial 
expansion, and then evaluating the polynomial, we need 


(n® — 2n + 1) + kn M-D, 

n(n — 1)? + kn A-S. 
To obtain k values of D(A) from the determinant will require 
(k/3)(n* + 2n — 3) M-D, 

(kn/6)(n + 1)(2n + 1) A-S. 


| A comparison of these results shows that it will be quicker to obtain the polynomial 
expansion first if we need more than three values of D(A). 
If we use iterative methods to solve (1), we form the sequence of matrices 


AX, A®X,---, (2) 
where X is an arbitrary column matrix,* 
X = {X;,+--, Xe}. 
To form each member of the sequence (2) requires m* multiplications and divisions 
and n(n—1) additions and subtractions; hence for k iterations we need 
kn? M-D, 
kn(n — 1) A-S. 
If we have the polynomial form, each iteration requires only ” multiplications 


and divisions and »—1 additions and subtractions. Adding these to the operations 
required to expand (1) into polynomial form by the modified Danielewsky method 


a (§2c) we need 
(n? — 2n + 1) + kn M-D, 


be n(n — 1)? + k(n — 1) A-S. 


Hence if k2n+1 the expansion to polynomial form represents a net saving, still 
using the powerful iterative method [3, 4, 16, 17, 18]. 


II. EXPANSION OF A DETERMINANTAL EQUATION INTO POLYNOMIAL FORM 


2. Methods applicable to the case | A — JA| =0. Direct expansion is tedious except 
for the very lowest orders, although sometimes it is desirable because all the elements 
need not be given numerically. Purely numerical methods, such as the method of un- 
determined coefficients or the use of an interpolation formula, will be described further 


* In this paper, a row of quantities enclosed in braces will denote a column matrix. 
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on, as they are applicable to the most general case. Of the various methods particu- 
larly applicable to the present case (Eq. 1), five will be discussed. These will be taken 
up in chronological order of discovery. 

a) The method of Leverrier [7, 12, 19]. Until recently, Leverrier’s method was 
probably the best general method for obtaining the polynomial expansion of the char- 
acteristic equation of the matrix A. Let the characteristic values of this matrix, i.e., 
the roots of the equation =0, be Au, Ae, , Then from well known rela- 
tions between the coefficients of a polynomial equation and its roots [20], we have 

|A — = (— + (— + 
+ + (3) 
But direct expansion of the determinant gives 
| A — = (= + (= 1)" Maar + ang) 
whence we have the relation 
+ doe + +++ + = Ar (4) 
If we consider the kth power of the matrix A, the characteristic values of the matrix 
A* will be Xj, M3, - , [21]. From this we obtain 
nn n— k k n— 
(k) (k) 


= (— + (— + + Onn + ore, 


which yields the relation 


(k) (k) 


where the terms a are taken from the principal diagonal of A*. 
If we write our original equation as 


|A — = (— + + + + + pa] = 0, 
we have, from (3), (4) and (5), 
pi = — (Art +++ +n) = — (11 + + + Gan) = — Sr, 
2 


pi = — sipi = 


or 


pi = — Si, = — + 


which are of a set of simultaneous linear equations for the coefficients p, (k=1, 


+, 
Leverrier’s method can be summarized as follows: We form the powers of the 
matrix A, i.e., A* (k=1, 2, - +--+, m), and add the diagonal terms of each matrix to 


obtain 


| 

| 

a 

| 
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(k) (k) (k) 


Sk = Qyy + + Onn. 
We then set up the set of m simultaneous linear equations 
pe = — (1/k)(Sipe_1 + +++ + se) (k=1,2,---,m), (6) 
which can be solved for p; (R=1, 2, - - - , m). We can then write the equation 
D(x) =| A — D| = (— + + + + pp). 


Leverrier’s method is rather tedious because of the labor of forming the powers 
of the matrix A. Each element of a product matrix will require m multiplications and 
n—1 additions and subtractions, and each matrix will have n* terms, except that one 
need form only the » diagonal terms for the last matrix. The fastest way to solve the 
simultaneous equations is to solve them successively, i.e., to use the solution of the 
first to solve the second, the solutions of the first two to solve the third, and so on. 
This will require 4(m?-++n—2) multiplications and divisions and $n(m—1) additions 
and subtractions. Consequently, Leverrier’s method will require, in the general case, 


3(m — 1)(2n* — 2n? + n + 2) M-D, (7) 
an(n — 1)(2n? — 4n + 3) A-S. 
Example. Let us consider the matrix 
6 —3 4 1 
4 2 + 0 
A= 
4 —2 3 1 
4 2 3 1 


The sums of the elements of the principal diagonals of A, A*, A*, A‘ are found to be 
$1=12, s2=56, 53= 288, ss=1504. Thus Eqs. (6) take the form 


pi = — 12, 
p2 = — (1/2)(12p; + 56) = 44, 

ps = — (1/3)(12p2 + 5661 + 288) = — 48, 

ps = — (1/4)(12p3 + 56p2 + 288p1 + 1504) = 16, 


and we have 
| A — = — + — 48 + 16. 


b) The method of Krylov [7]. Even as originally formulated by Krylov, this method 
represents a considerable saving of effort over the method of Leverrier when the order 
of the determinant is greater than four. As modified by Fraser, Duncan and Collar 
[22], the saving is even greater. The modified form of this method is as follows. 

The Cayley-Hamilton theorem [23] states that a square matrix satisfies its own 
characteristic equation when interpreted as a matrix equation, i.e., if 


+ + pre? +--- +p, = 0 
is the characteristic equation of the matrix A, then 


A*™ + + poA™?+---+ = 0. (8) 
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If we postmultiply (8) by an arbitrary column matrix 
C(0) { cro, Cro}, 
and define the sequence 


C(1) AC(0) { cu, Cay ca}, 
C(2) = AC(1) = A*C(O) = Gna}, 


we obtain the matrix equation 
C(0)pn + + + C(m — 2)p2 + C(n — 1)pi = — C(n). 


This is equivalent to the set of m simultaneous linear equations in the » unknowns 


Pi, pe, » Pn; 


n—1 


CikPn—k = — Cin (¢ = 1,2,---,m). (9) 
k=0 
Solving this set of equations for the p,’s, we can readily write down the polynomial 
equation. 

The formation of each of the column matrices C(k) requires n? multiplications and 
n(n —1) additions and subtractions. The solution of Eqs. (9) by Aitken’s method [28] 
requires (1/2)n?(n+3) multiplications and divisions and (n/2)(m?—1) additions and 
subtractions. Consequently, we require 


(3/2)n*(n +1) M-D, (10) 
(n/2)(n — 1)(3n + 1) A-S. 
Krylov’s original method can be shown to require 
(1/3)(n* + 4n? + 2n? — n — 3) M-D, (11) 
(1/6)n(n — 1)(2n? + 7n —1) A-S. 


Example. Let us again consider the matrix of the previous example. We have 
C(0) = {1,1,1,1}, C(1) = {8,10,6,10}, C(2) = {52, 76, 40, 80}, 
C(3) = {324, 520, 256, 560}, C(4) = {1968, 3360, 1584, 3664}. 


Thus (9) becomes 
Pat 8p3 + S2p2 + 324p; + 1968 = 0, 


pa t+ 10p3 + + 520f1 + 3360 = 0, 
Pat Ops + 40p2 + + 1584 = 0, 


whence 
ps = 16, ps = — 48, po = 44, fi = — 12. 


The characteristic equation of A is then 


16 — 48 + 44d? — 124° + Af = O. 


| 

q 
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(k) (k) (k) 


Sp = + dog 
We then set up the set of » simultaneous linear equations 
pe = — (1/k)(Sipe-1 + +--+ +Serpit se) (k= 1,2,---,m), (6) 
which can be solved for p; (k=1, 2, - - - , 2). We can then write the equation 
=| A — = (— + pir! + por? + + pp). 


Leverrier’s method is rather tedious because of the labor of forming the powers 
of the matrix A. Each element of a product matrix will require m multiplications and 
n—1 additions and subtractions, and each matrix will have n* terms, except that one 
need form only the ” diagonal terms for the last matrix. The fastest way to solve the 
simultaneous equations is to solve them successively, i.e., to use the solution of the 
first to solve the second, the solutions of the first two to solve the third, and so on. 
This will require 4(m?+-n—2) multiplications and divisions and }n(n—1) additions 
and subtractions. Consequently, Leverrier’s method will require, in the general case, 


4(m — 1)(2n* — 2n? + n + 2) M-D, (7) 
3n(n — 1)(2n? — 4n + 3) A-S. 
Example. Let us consider the matrix 
6 -—3 + 1 
4 2 4 0 
A= 
4—2 3 1 
4 2 3 1 


The sums of the elements of the principal diagonals of A, A?, A*, A‘ are found to be 
$1=12, se =56, s3= 288, ss=1504. Thus Eqs. (6) take the form 


fi = — 12, 

p2 = — (1/2)(12p1 + 56) = 44, 

ps = — (1/3)(12p2 + 56p1 + 288) = — 48, 

ps = — (1/4)(12p3 + 56p2 + 288, + 1504) = 16, 


and we have 
| A — = — + 44? — 48 4+ 16. 


b) The method of Krylov [7]. Even as originally formulated by Krylov, this method 
represents a considerable saving of effort over the method of Leverrier when the order 
of the determinant is greater than four. As modified by Fraser, Duncan and Collar 
[22], the saving is even greater. The modified form of this method is as follows. 

The Cayley-Hamilton theorem [23] states that a square matrix satisfies its own 
characteristic equation when interpreted as a matrix equation, i.e., if 


+ + +p, = 0 
is the characteristic equation of the matrix A, then 


A” + + +--+ + pal = 0. (8) 
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If we postmultiply (8) by an arbitrary column matrix 
C(0) { Cro}, 
and define the sequence 


C(1) = AC(0) = { ** cu}, 
C(2) = AC(1) = A°C(O) = {c12, Cao}, 


we obtain the matrix equation 
C(0)pn + C(1)pn-1 + + C(m — 2)p2 + C(n — = — C(n). 
This is equivalent to the set of simultaneous linear equations in the unknowns 


Pi, Po, » 


n—1 
> CikPn—k = — Cin (¢ = 1,2,--+,m). (9) 
k=0 

Solving this set of equations for the p,’s, we can readily write down the polynomial 


equation. 
The formation of each of the column matrices C(k) requires m? multiplications and 


n(n —1) additions and subtractions. The solution of Eqs. (9) by Aitken’s method [28] 
requires (1/2)n?(m+3) multiplications and divisions and (n/2)(n?—1) additions and 
subtractions. Consequently, we require 


(3/2)n?(n +1) M-D, (10) 
(n/2)(n — 1)(3n + 1) A-S. 
Krylov’s original method can be shown to require 
(1/3)(n* + 4n? + 2n? — n — 3) M-D, (11) 
(1/6)n(m — 1)(2n? + 7n —1) A-S. 


Example. Let us again consider the matrix of the previous example. We have 
C(O) = {1,1,1,1}, C(1) = {8,10,6,10}, C(2) = {52, 76, 40, 80}, 
C(3) = {324, 520, 256, 560}, C(4) = {1968, 3360, 1584, 3664}. 


Thus (9) becomes 
pat 8p3 + 52p2 + 324p1 + 1968 = 0, 


pa + 10f3 + 76p2 + 520f; + 3360 = 0, 
pat Ops + 40p2 + 256p; + 1584 = 0, 


ps + 10p3 + 80p2 + 560f; + 3664 = 0, 


whence 
Ps = 16, ps ae 48, pe = 44, pi soe 12. 


The characteristic equation of A is then 


16 — 48 + 44d? — 12d + = O. 


| 
| 
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c) The method of Danielewsky [8]. The essence of Danielewsky’s method is the 
transformation of the expression D(A) =|A—JA| to the Frobenius standard form 


ps Pn 
1 — 0 0 

D()= | 0 1 O | --- —p,). 
0 0 — 


This gives the polynomial expansion directly. 

Danielewsky starts with the (n—1)th element of the mth row (a constant term), 
reduces it to unity, and then uses this to eliminate the constant terms from the other 
elements of the mth row. This process introduces extraneous terms in A in the (n—1)th 
row, which can then be removed by multiplying the other rows by appropriate con- 
stants and adding to the (n—1)th row. A similar procedure is then followed with the 
(n—2)th element of the (x —1)th row, and the reduction is continued until the stand- 
ard form is reached. 

This process of elimination can readily be carried out by matrix multiplications. 
Let us consider a matrix of order 6 which has already had two rows reduced. It is 
then of the form 


he Ci2 C13 Cis 
Co1 Coz Cog Cog Cop C26 
C31 C32 C33 C34 C35 C36 
C= 
Car (Can Cap Cag Can Cae 
0 O 1 
Loe ¢ 1 0 


and the determinant has been reduced to the form 


Di) =|C— 


We now postmultiply the matrix C by the matrix E, 


i 1 0 0 0 0 0 7 
0 1 0 0 0 0 

—Ca2/ Cas 1/c4s —Cas/Caz 
0 0 0 i 0 0 
0 0 0 0 1 0 

| 0 0 0 0 0 a 


to obtain 


f 
= 
| 
| 
| 4 . 
4 
a | 
| 
| 
4 
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Cir Ci2 Cig Cia Crp C16 
Cos Cos Cop C26 
Cc’ C31 C32 Ces C35 C36 
£- 
reo. 


After this transformation we have 
=|C’ — - 


If this expression is now premultiplied by E-!, we return to a form which is equal ‘to 
our original expression, but one step closer to the standard form, 


=| — =|C” — |. 


Fortunately, E~! can be written down directly, 


Car Can Cap Cag Cas 


This particular premultiplication changes only the third row of C’, hence the determi- 


nant has been transformed to the form 


Cu — C12 C13 Cu C15 C16 

C21 — X C23 Co C5 C26 

- 6 

D(r) = 31 32 33 34 35 3 
0 0 1 —rX 0 0 

0 0 0 1 —r 0 

0 0 0 0 1 —X 


A continuation of this process will eventually yield the normal form. 

If the matrix method is followed strictly, it involves an undue amount of writing. 
With only a slight increase in the number of operations, it can be abridged to give 
greater ease of calculation with a calculating machine, and to permit checking at every 
stage of the computation. A numerical example will best illustrate the method and 


the check. 


Example. For the matrix of the two previous examples, the scheme of calculations 
will run as follows: 


q 
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(1) 
(2) 
(3) 
(4) 


(4’) 


(S) 
(6) 
(7) 
(8) 


(9) 
(9’) 


(10) 
(11) 
(12) 
(13) 


(14) 

(15) 
(16) 
(17) 
(18) 


(19) 


(14’) 


48 

11.3333 
—44 

45 .3333 


Hence 
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6 —3 4 1 8 

+ 2 4 0 10 

4 —2 3 1 6 

4 2 3 1 10 

(1.3333 0.6667 1 0.3333 3.3333) 

0.6667 —5 .6667 1.3333 —0.3333 —4 —5 .3333 
— 1.3333 —0.6667 1.3333 —1.3333 —2 —3 .3333 

0 —4 1 0 —3 —4 

0 0 1 0 1 

0 — 36 12 —4 —28 

(0 1 —0 .3333 0.1111 0.7778) 

0.6667 0.1574 —0.5556 0.2963 0.5648 0.4074 
— 1.3333 0.01852 1.1111 — 1.2593 — 1.4630 —1.4815 

0 1 0 0 1 

0 0 1 0 1 

48 11.3333 —44 45 .3333 60 .6667 

(1 0.2361 —0.9167 0.9444 1.2638) 

0.01389 0 0.05556 —0.3333 —0.2639 —0.2778 

1 0 0 0 1 

0 1 0 0 1 

0 0 1 0 1 

12 —44 48 —16 0 


D(A) = — + 44? — 480 + 16. 


The explanations of this scheme are as follows. We first postmultiply the matrix 


1 

0 
—4/3 

0 


0 0 

1 0 
—2/3 1/3 

0 0 


0 
0 


1 


whose elements are given in lines 1, 2, 3 and 4, by the matrix 


—1/3 


This is accomplished by dividing the elements of row 4 by the element in the third 


column, 3, yielding row 4’. The second-order minors of the unit element in the third 
- column of row 4’ are then formed with rows 1, 2, and 3, the unit element always 
; taking the leading position. This is easily done by writing row 4’ on a card, and form- 
ing the cross products with the various rows. These minors are entered in rows 5, 6 
and 7, under the column corresponding to the elements with which the cross product 


is formed. The elements in column 3 are formed by dividing the corresponding ele- 


be: ments of column 3, rows 1, 2 and 3 by the italicized element in row 4. Row 8 is im- 
f ; mediately written as shown. Thus, the element — 5.6667 in row 5, column 2, comes 
from 1(—3)— (0.6667)(4), and the element 1.333 in row 5 column 3 comes from divid- 
ing the element 4 in row 1 column 3 by 3. 


| 
|_| 
| 
= 
—36 
3 12 
—4 
I 
x 
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We next premultiply the matrix with the elements given in rows 5, 6, 7 and 8, by 
the matrix 


1000 
0 1 0 

4231 
0001 


This is accomplished by writing the elements of row 4 in a column on a card (here 
shown written to the left of rows 5, 6, 7 and 8), and forming the sum of the products 
of these numbers with the elements of the columns of the rows 5, 6, 7, and 8. This 
yields row 9, which is the transformed form of row 7 after the matrix multiplication. 
Since the rest of the matrix is unchanged, it is unnecessary to rewrite it. 

The whole process is now repeated starting with row 9, dividing each element of 
that row by the italicized element (—36) to obtain row 9’. This is written on a card 
and used to form the cross products with rows 5 and 6, giving the elements in the 
first, third and fourth columns of rows 10 and 11. The elements in the second column 
are obtained by dividing the corresponding elements of rows 5 and 6 by —36. Rows 12 
and 13 can be written down immediately as shown. The process is continued until row 
19 is reached. At this stage it is unnecessary to rewrite the entire matrix, since the 
desired coefficients appear in only the first row, i.e., row 19. The polynomial can row 
be written down as shown. 

The columns labelled }> and >-’ are used for checking the work. The elements 
in the column labelled >> are obtained by summing the elements of the first four col- 
umns of that row, while those in pm come from only three columns, omitting that 
column which contains the element used as the pivot for the previous set of cross 
multiplications. The cross products formed with the >> columns should equal the ele- 
ments of the } ny column at the next stage of the transformation, e.g., the element 
— 3.333 in row 6, column >’ comes either from adding the elements (—1.3333) 
+(—0.6667)+(—1.3333) of row 6, or from the cross product (1)(10) — (4)(3.3333). 
Since these give equal results, the computation of row 6 is probably correct. This 
check is not applicable to the row just reduced, so there is no point in forming >,’ 
for that row, or the rows already in standard form. The accuracy of row 9, and similar 
rows, is checked by forming the sum of the products of the column (4, 2, 3, 1) and the 
elements of column >... Since this product-sum is equal to the sum of the elements of 
row 9, the accuracy of that row is checked. Compensating errors can occur, so the 
check is not absolute, but it is a great help in avoiding an accumulation of errors. 

We must next consider the exceptional case in which a zero appears for the ele- 
ment with which we expect to divide in making the next reduction, i.e., the element 
one place to the left of the diagonal. The following two cases arise: (1) There is at 
least one element in the row under consideration which does not have a vanishing 
constant term. (2) All of the constant terms in the row under consideration vanish. 

Case (1) can be decomposed into two subcases, according as the non-vanishing 
element is (a) to the left of the diagonal, (b) to the right of the diagonal. In subcase 
(a), we add the elements of the column containing the non-vanishing element to the 
column in which we wish to introduce a constant term. (This technique can also be 
used if the element immediately to the left of the diagonal has a fairly large tabular 
error, and another term farther to the left is more certain.) This will not only intro- 


| 3 4 
4 
| 
| 
| 4 
J 
| 
| 
| 
A 
= 
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duce the desired constant term, but in some row it will introduce an unwanted term 
in \ off the diagonal. This can be removed, however, by subtracting the appropriate 
row from the row containing the extraneous \. The reduction can then go ahead as 
usual. In subcase (b), the determinant is immediately factorable into the product of 
two determinants, one of which is already in standard form. 

These subcases can best be illustrated by examples. In subcase (a), let us suppose 
that after two reductions we reach the form 


3 —2 5 3 
1 2—-\A -1 4 1 
2 0 4—-ir 6 |. 
0 0 1 — 0 
0 0 0 1 ae 


If we add column 1 to column 2, we obtain 


7—-\’ —2 5 3 
1 + 1 
2 2 4—-r -1 6 
0 0 1 — 0 
0 0 0 1 —\ 


This has an extraneous A in row 1, which we can eliminate by subtracting row 2 
from row 1, to obtain 


4 -1 1 2 
1 4 1 
2 2 4—i -1 6 
0 0 1 —A 0 
0 0 0 1 —> 


This is now in a form capable of treatment by the general method. In subcase (b), 
we might reach the form 


4— 3 —2 5 3 
1 + 1 
0 0 -1 6 |, 
0 0 1 —X 0 
0 0 0 1 —A 
which can be factored into the product 

-1 6 
3 

1 — 
1 2— 

0 1 


q 

| 
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The determinant on the right is already in the Frobenius standard form, while that 
on the left can be expanded immediately, although in the general case it would have 
to be reduced further by the general method. 

In case (2), the vanishing of all constant terms in a given row indicates that ) is a 
factor of D(A). If the rth row from the bottom has vanishing constant terms, it means 
that A’ is a factor of D(A). The determinant for the lower degree polynomial which 
we have yet to determine can readily be constructed from the elements above the 
vanishing row. As an example, let us consider 


3 —2 5 3 
1 2-A -1 4 1 
0 0 0 0}. 
0 0 1 — 0 
0 0 0 1 —A 
This is equal to 
4- 3 | (— 
1 
In its original form, Danielewsky’s method requires 
(n? — 2)(m — 1) M-D, (12) 
n(n — 1)? A-S, 
and in the modified form given in detail above, it requires 
(n — 1)(n? + n — 1) M-D, (13) 
n(n — 1)? A-S. 


In spite of the extra operations required, the modified form is to be preferred, because 
it is better adapted to routine computation with a calculating machine, and because 
it can be checked at each stage of the computation. 

d) Reiers¢l's method. Reiers¢l [14] bases his method of obtaining the coefficients 
of the determinantal equation 


D(d) =| 4 — DA] = (— — — — — pa) 


on the fact that the coefficients p; can be calculated as (—1)*+! times the sum of all 
k-rowed principal minors of the matrix A. The method is powerful for low values of n, 
but for large m the labor is considerable. 

Reiersgl uses a method for computing the principal minors of A based on the 
same pivotal method used by Aitken in various numerical processes dealing with de- 
terminants [28]. In the process of evaluating a determinant by Chid’s method [15], 
simple quotients of various minors are obtained, and the method is easily extended to 
give all of the principal minors of the matrix. 

This form of calculation is used here since it is uniform with most of the other 
methods described in this paper, and requires essentially the same number of opera- 
tions as the application of Reiers¢gl’s recursion formulae. In fact, it is merely a schema- 
tization of his method. 


| 
j 
it 
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In thie process of evaluating the determinant |A| by Chid’s method [15], we use 
the identity 
Gig 413° ** Ain 
Go2 423° Aan 
G32 «33° * * 


| 
|A| = G33°** Gan | = C11 
, 
Qn2 Gn3°** Ann 
| Qn1 Gna Gnz*** Ann 
where aj, (i, 7=2, 3, ---,m). If we define our principal minors (mi- 
nors obtained by striking out the same rows as columns) as B,,.... (r<s<--- <w), 


where 7, s, - - + , w are the indices of the rows (and columns) used to form the minor, 
then we have 


, 
Bi, = 2114s. 


Carrying the reduction one stage farther, we obtain 


ut 


G33 °° Gin 
4? 
44° * Dan 
| A | = @11022 
Ann 


where, as before aj/ = aj, Then 
47 
= 411022411, 


and so on through the reduction. This gives all the second order principal minors of 
the form B;,, the third order minors of form Bi;, the fourth order minors of the form 
Bysu, and so on, including the value of the determinant itself . 

For the minors in which the first two subscripts are 1 and 3, we start with the 
determinant 


a33 * G3n 
G43 * Aan 
’ 
Qn3°*** Onn 


and carry through the pivotal reduction as before. For minors of the form Bz... we 
start with 


Gan 
G32 433° Gin 
Ann 


In this way all of the principal minors can be built up. Since we always start with a 
term on the principal diagonal, the algebraic sign presents no problem. 
The number of operations required by this method is: 


| 

a 

F 
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5-2" — (n?+ 4n+ 5 M-D, ere 
(n? + 4n + 5) (14) ee 
4-2" — (n? + 3n + 4) A-S. 
Example. Let us again consider the expression 
DO) 4 4 0 
3-r 1 4 
4 2 3 
The scheme of calculations will run as follows: 4 
an = 6 2/3 1/6 B, = 6 
4 4 0 B, = 2 
4 —2 3 1 B; = 3 x 
4 2 3 1 
0 1/3 1/3 By = 6-1/3 = 2 q 
a3; = 1/3 1 1 Bus = 6-4-1/3 = 8 2 
—1 1 Byy = 6-4-1 = 24 
2 1 Bie = 6-4-(1/3)-2 = 16 
aq = 6 
ai; = 1/3 1 1 
1/3 1/3 
0 Bim = 6-(1/3)-0 = 0. 


This completes all terms with 1 as the first index. Starting with the third order de- 
terminant we obtain by striking out the first row and column, 


de = 2 1 2 0 
—2 3 1 
2 3 1 
bis = 7 1/7 Bos = 2-7 = 14 
oun 1 Bu = 2-1 = 2 
b!” = 8/7 1 Boy = 2-7-8/7 = 16 


\ 

3 1 

Pu = = 0, 

3 1 

4 
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Hence we have 
pi = (— 1)9(Bi + Be + Bs + By) = 12, 


bo = (— 1)8(Bire + Bis + Bis + Bos + Bas + Bus) = — 44, 
pe = (— 1)*(Bizz + Bisa + Bisa + Boss) = 48, 
pa = (— 1)5-16 = — 16, 

D(d) = 4 — + 44d? — 48 + 16. 


e) Samuelson’s method. Samuelson [13] has devised one of the fastest methods yet 
developed. His method requires a few more operations than Danielewsky’s, but the 
routine involved is extremely simple. 

To get the polynomial expansion of 


=| A — = (= — pid! — — — 
we consider the differential equation 


(n) 


D(d/dt)x,() = (— — pint” — pom” — — = 0, (15) 


where the superscripts in brackets denote derivatives with respect to ¢. Equation (15) 
can be written as a set of m simultaneous first order differential equations 


where 
1 0 


0 1 O---0 O 


is the “companion matrix” to the polynomial in question. 
Actually we need a scheme to go from a system in many variables to a high order 
system in one variable. Samuelson accomplishes this in the following manner. 
Let us consider the system 
Ax(t) = x'(d), (16) 
where x(f) is the column matrix 


x(t) = {a(t), 


Equation (16) gives us equations in the 2m variables x1, x2, + ,%n,%1,%2, °° 
There are insufficient equations to eliminate all of the variables except those carrying 
the subscript 1. However we can differentiate (16) (~—1) times with respect to ¢, 
obtaining the nm? equations 

Ax (2) x™ (2), 


= an) 


Ax(t) = 


We now have linear equations in the m?+ variables (x1, x2, 


4 
4 
{ 
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x). We can use all but one equation to eliminate the 


n?—1 variables not involving the subscript 1, and substitute in the remaining equa- 
tion to get the desired high order equation in x; and its derivatives. 
Let us consider 


a1 R 
Gon | = 
A= 
Int Gn2*** Ann 


If we transfer the variables with subscript 1 to the right of (17), we can rearrange 
and rewrite it in the form 


0 M:-- 0 0 0 -S --- 0 


0 0 —I M; 0 0 0 -S§ 
0 0 O R 0 0 
. (18) 
0 0 O--- R 0 0 0 0 --- O 
0 0 R 0 0 0 1 —~@Ges 0 


0 R O 0 0 1 O 


The elimination of the »?—1 unwanted variables from the first »?—1 equations and 
the subsequent substitution in the remaining equation can be performed by pivotal 
reduction [28], always using elements of the matrix on the left of (18) until a single 
row remains on the right. 

Reduction down to the first row containing R can be made in the general form, 
yielding 


0 0 0 0---0 0 1 
RM 0 0 0 0---0 1 -—ay, —RS 

RM? 0 0 0 O---1 —RS —RMS }. (19) 


In practice this is the point at which to start the reduction. 
To set up the matrix (19) requires 


n(n — 1)? M-D, 

n(n — 1)(n — 2) A-S. 
Pivotal reduction of (19) will require 

4n? — 13m + 12 M-D, 


4n? — 13n + 12 A-S, 


4 
‘ 
. 
ie 
? 
f 
- 
« 
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making totals for the method of 
n* + 2n? — 12n + 12 M-D, 
n?— 11n+ 12 A-S. 


Samuelson uses a method due to Crout for the reduction of his equations. Crout’s 
method involves forming exactly the same products and sums as are formed in 
Aitken’s method used above, although Crout’s formulation involves somewhat less 
writing than the above method, but also requires keeping in mind somewhat more 
complicated formulae. For the average engineer or physicist, ease is fully as important 
as speed. 

Example. We again consider the matrix 


6 |-3 4 
— an R 
2 4 0 j=| 
4 |-2 3 1 

2 3 


We find that RM=[—12, 3, 5], RM?=[-—20, —24, 8], RM*=[24, —128, —16], 
RS= [8], RMS=[—16], RM?S=[—144]. The matrix to be reduced then becomes 


—3 4 1 | 0 0 0 1 —6 
—12 3 5 | 0 0 i —6 —§ 
—20 —24 8 | 0 i -6 -8 16 

24 -6 -8 16 144 


The reduction then proceeds as follows: > 


(20) 


—3 -—-—/0 0 0 2 1 


| 
1 1 
| 3 5| 0 0 1 —6 —8 | —17 
—20 —24 8 | 0 1 —6 —8 16 | —33 
24 —128 1 —6 —8 16 144| 27 
1 1 10 16| 5 
| —13 i .-—|/0 ima 
oe 13 13 13 13} 13 
1 —6 56 | —13 
3 3 
| —-9%  -8/1 —6 —8 24 96 3 
1300 | o 5018 _ 12324 3224 | 3289 
a4 507 | 1300 1300 1300 1300 | 1300 
| 200 | 200 1272 288 | 519 
| 13 13 13 | 13 
1 44 —48 16| 1 
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The pivotal element for each succeeding reduction is made equal to unity by dividing 
that row by the value of that element. The column marked )_ is used as a check. For 
any stage of the reduction, the cross products are formed using the >, column as if it 
were any other column, and the values entered as usual. These values should equal 
the sum of the elements in the row in which they appear. The check is not absolute, 


but it is very useful. 
3. Methods applicable to the case |A — By| =0. -a) The method of reciprocation. 


The equation 
|A — Bar| =0 


has the same roots as the equation [24] 
| B14 — = 0, (21) 


where B-' is the reciprocal of B, provided only that B is not singular. The matrix 
product B-!A can readily be formed by Aitken’s method [28], and the determinant 
(21), which contains \’s only along the principal diagonal, can then be expanded by 
one of the methods of §2. 

_ The formation of the product B-'A requires 


(n?/2)(3n — 1) M-D, 
(n/2)(n — 1)(3n — 1) A-S. 
Using the modified Danielewsky method to obtain the polynomial form, we shall need 
all told 
(1/2)(m + 1)(5n? — 6n + 2) M-D, (22) 
(n/2)(m — 1)(Sn — 3) A-S. 


Example. Let us consider the determinant 


—9+2r —8+3X —7+d —74+2A 
15—3 16—S5A 13-2 


D(\) = 
—8+r —8+2\ —7+2\ —8+3A 
24-5 19-3 22—6d 
We have 
15 16 13 15 —3 -5 -2 -4 
—8 —8 -7 -8 et 
| 23 24 19 22 | -3 
0-2 37 0 3 
5 0 1-2 5 —2 
4 8 4—-r 7 


_ 

? = 
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1 2 -—-8 0-2 3 —9 —8 -—7 —7 
—5 -—4 1 1§ 16 13 15 
B(B-1A) = = =A, 
1 2 4 8 4 —8 —8 —7 
—3 —5 -—3 -6 —6 —8 —5 -—7 23 24 19 22 
Thus 


= + + 33A2 + BA + 8. 


The calculation of B-! can be carried through conveniently by means of Aitken’s 
method [28] of obtaining the reciprocal of a matrix. The product B(B-'!A) can be 
formed as a check. 

b) The Danielewsky-Masuyama method. An extension of Danielewsky’s method 
of transforming a \ determinant to the Frobenius standard form has been made by 
Masuyama [6] for the case | A —Bd| =0. This method requires 


(1/24)n(m — 1)(7n? + 13n + 66) M-D, 


(23) 
(1/24)(m — 1)(7n® + Sn? + 58n — 48) A-S. 


For low orders (through the fifth) this represents a small saving in the number of 
operations over the method given in §3a, but it is not as well adapted to machine 
computation. It is applicable, however, to the case in which the determinant of the 
matrix B vanishes, but in this case the method of §4b, using Newton's interpolation 
formula, is to be preferred. For these reasons, we shall not consider the method in 
more detail here. 

4. Methods applicable to the case | AA*+A,A""+ --- +A,| =0. -a) Transfor- 
mation to the form | A —l| =0. It is possible to transform an mth order determinant, 
the terms of which are polynomials at most of degree m in X, into a determinant of 
order mn with terms linear in A, provided that the matrix of the coefficients of \” is 
not singular. This can be done in more than one way, but the following seems most 
convenient [25]. 

If the determinant we wish to transform is 


D(d) = | + +--+ + + A, | = 0, (24) 
we consider the related set of simultaneous linear differential equations 
(ApD* + = 0 (D = d/dt), (25) 
where x is the column matrix 
x= { x, ° °° » Xn}, 
and 
| Ao| ¥ 0. (26) 


Because of the condition (26) the reciprocal matrix A! exists. Consequently we can 
premultiply Eq. (25) by Ag}, obtaining 


(ImD* + + --- + + = 0, (27) 


where J,, is a square matrix of order m with units on the principal diagonal and zeroes 
everywhere else. We can write (27) in the form 


ImD + Ag 4+ + + + = 0, (28) 


= 
4 
; 
j 
: 3 
| 
| 
F 
4 
— 
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where the superscripts in brackets denote derivatives with respect to ¢. Equation (28) 
represents a set of m simultaneous differential equations. There are n of the column 
matrices x“. Thus to make the set equivalent to Eq. (27) we need the m(n—1) 
additional equations of the form 


a) = = = = Dg = 1,2,---,m— 1). (29) 


If in Eqs. (28), (29) we now set 


we get the following set of mn homogeneous simultaneous linear equations in the g’s: 
— Imd)g**—? + Ag + +++ + + = 0, 
= = 0, 


— — Ingdr = 0. 


The condition of compatability for this set of equations is the vanishing of the com- 


pound determinant 
Ag'A, Ap "Ana 
TI. 0 0 
0 —In 0 
0 0 --+ —In — Imd 


Since the \’s appear only along the principal diagonal, this can be expanded into poly- 
nomial form by any of the methods of §2. 

This method involves the calculation of the reciprocal of Ao and the formation 
of m matrix products of the form A~!B. This can be done by Aitken’s method [28]. 
The special form of the expression makes the expansion to polynomial form by 
Danielewsky’s method considerably faster than in the general case. 

In the quadratic case, which is the most general usually met with in physical prob- 
lems, the transformation to diagonal form requires 


3m* M-D, 
m(m — 1)(3m — 1) A-S, 
while the reduction to polynomial form requires 
2(3m — 1)(m* — 1) M-D, 
m(6m* — 9m + 5) + 2 A-S. 
The total number of operations in the quadratic case will be 
9m? — 2m? — 6m + 2 M-D, 
9m* — 13m? + 6m : 2 A-S. (30) 


Example. Let us consider the determinant 


é 
| 

= 
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M4 O54 R43 OF 4 3 
+2 + 4d? + 3d 
We have 
ri -1 45 5 2 15 
0 Awl 1 #44 
2-1 O- 3 4 3-7 
10 —11 57 30 24 ~—59- 
if 
46 46 
2 7 1 
i$ 13 50 507 
if 
-65 S54 87], —41 -—94 —48 
46 46 7 


Consequently we can write our determinantal equation in the diagonal form | A — JA], 
where A is the square matrix 


j 30/46 24/46 —59/46 15/46 —16/46 —13/46 13/46 50/46 50/467 
—176/46 —12/46 —5/46 —65/46 54/46 87/46 —41/46 —94/46 —48/46 
6/46 14/46 25/46 3/46 6/46 25/46 21/46 10/46 10/46 
—1 0 0 0 0 0 0 
0 ~! 0 0 
0 0 ~!1 0 0 0 0 0 
0 0 0 -1 0 0 0 0 0 
0 0 0 0 | 0 0 0 0 
a 0 0 0 0 0 =~] 0 0 o J 


b) Interpolation-formula method. The mth order determinant (24) will, on expan- 
sion, give a polynomial of degree not greater than r=mn in X. Such a polynomial 
will contain (r +1) numerical coefficients of the various powers of \, and if we evaluate 
D(A) for (r +1) numerical values of X, we will have sufficient data to determine the 
coefficients. This can be done either by using an interpolation formula or by solving 
a set of simultaneous linear equations for the coefficients. 

If the successive values chosen for \ differ by a constant difference, the Gregory- 
Newton interpolation formula is well adapted to obtaining the polynomial expansion. 
This requires the formation of a difference table. However, the difference table is 


= 
» 
’ 
4 
as 
+ 
4 
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easy to compute, and in addition it indicates the order of polynomial to expect, which 
is particularly useful when the matrix of the coefficients of X" is singular (| Ao| =0). 


The Gregory-Newton interpolation formula is usually written in the form [26] 


= D(a) + xAD(a) + (0) 


r! 


where \=a+<xh, and the difference functions A’D(a) are formed as follows: 


D(a) 
AD(a) 
D(a + h) A*D(a) 
AD(a + h) A*D(a) 
D(a + 2h) A’D(a + h) 
AD(a + 2h) 
D(a + 3h) 
D(a + rh) 


where AD(a) = D(a+h)—D(a), A?D(a) =AD(a+h) —AD(a), and so on. The kth differ- 
ences of a kth degree polynomial will be constant, so this gives us immediately the 
degree of the polynomial. It will always pay to calculate one or two extra values of 
the determinant to check the constancy of the last differences as a check on one’s 
work. 

Equation (31) does not yield the polynomial form directly, since each term is a 
polynomial in \. The transformation to descending powers in \ can be made once for 
all, so that the polynomial form can readily be calculated once a difference table is 
worked out [11]. We wish to put D(A) in the form 


D(d) = — = pot — a) + pd — (32) 
t=0 
This can be done by using the relation 
A*D(a)a,(s) 


Po = D(a), 


where the terms a,(s) are defined by the equation 


4 
| | 
4 
~ 
| 
2 
r 
+ 
- 
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r! 
s! = > a,(s)r*. 


For the purpose of calculation it is easier to use the recursion relation 
ax(s + 1) = apr(s) — saz(s). 

Table I gives numerical values of the function 


bi(s) = a,(s)/s! 


which, when used with the tabular differences and Eq. (33) will give the polynomial 
equation very quickly in the form (32). If a has not been chosen equal to zero, the 
equation can readily be rearranged in descending powers of X, instead of \—a, by 
synthetic division. 

TABLE I. dy(s) =a,(s)/s! 


| 2 | 3 | 4 5 

1 | 1.00000 00000 | 

2 |—.30000  00000/+.50000 00000] | 

3 |+.33333 33333—.50000 00000|-+.16666 66667) 

4 |—.25000  00000|+.45833 33333|—.25000  00000/+.04166 66667 

5 |+.20000 00000-41666 66667|+.29166  66667|—.08333 33333 

6 |—.16666  66667/+.38055 55556 —.31250  00000+.11805  $5556|—.020833 33333 

7 |4+.14285 71429-35000  00000/+.32222 33333|+.034722 22222 

8 |—.12500  00000|-+.32410 71429|-.32569  44444/+.16788  19444|—.048611 11111 

9 |+.11111 11111]~.30198 41270)+.32551 +.061863 42593 
10 |—.10000 00000/+.28289 68245|—.32316  46825|+.19942 68078 ~.074218 75000 

she 6 7 | 8 9 10 


+.0013888 88889 
—.0041666 66667|+.0001984 12698 
+.0079861 11111)/—.0006944 44444/-+.00002480 158730 
—.0125000 00000)+-.0015046 29630|—.00009920 634921/+-.00000275 573192 
+.0174363 42593|— 0026041 66667 |+-.00023974 867725|—.00001240 079365|+-.00000027 55731922 


ONAN WHS 


— 


If the successive values chosen for A do not differ by a constant quantity, the 
Newton interpolation formula cannot be used. The Lagrange interpolation formula 
[27] is available, but the method of undetermined coefficients will usually prove more 
satisfactory. 

Both the interpolation formula method and the method of undetermined coeffi- 
cients require, first of all, the evaluation of r+1 numerical determinants of the mth 
order. This will require 


(r + 1)(1/3)(m — 1)(m* + m + 3) M-D, 


(r + 1)(1/6)m(m — 1)(2m — 1) A-S. 


: 

is 

ay 
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Each numerical determinant will require the evaluation of m? polynomials of de- 
gree n, each evaluation requiring m multiplications and m additions and subtractions, 
if we use the method of synthetic division. If we assume that 0 is chosen as one value 
of \, so that we need evaluate the elements of only r determinants, the total number 
of operations to obtain the r+1 values of D(A) will be 


mr? + (r + 1)(1/3)(m — 1)(m? + m + 3) M-D, 


(34) 
mr? + (r + 1)(1/6)m(m — 1)(2m — 1) A-S. 


The interpolation formula method will require an additional (}r)(r+1) subtrac- 
tions to form the difference table, and ($r)(r +1) multiplications and (4$7r)(r —1) addi- 
tions and subtractions to form the coefficients, if the difference between successive 
assumed values of X is 1. If not, we shall need an additional 7 divisions by powers 
of h, and r—1 multiplications to form the powers of 4. For the case h=1, the total 
for the interpolation method (including the evaluation of the numerical determi- 
nants) is 


3(2m + 1)r? + (r/6)(2m* + 4m — 3) + (1/3)(m* + 2m — 3) M-D, 


35) 
(m + 1)r? + (r + 1)(m/6)(m — 1)(2m — 1) A-S. ( 
Example. Let us consider the determinant 
DA) 4—3A+A7+ AF 2— A+ 3A? — 
| A! A? Ae as | as 
D(0) 20 
— 36 
D(1) —16 —476 
—512 — 2628 
D(2) — 528 — 3104 — 4968 
— 3616 — 7596 — 2880 
D(3) —4144 — 10700 — 7848 0 
— 14316 — 15444 — 2880 
D(4) — 18460 — 26144 — 10728 0 
— 40460 — 26172 — 2880 
D(5) — 58920 — 52316 — 13608 
— 92776 — 39780 
D(6) — 151696 — 92096 
— 184872 
D(7) — 336568 


The accompanying difference table shows that the polynomial will be of the fifth 
degree instead of the sixth. An investigation of the determinant of the coefficients of 
\* shows that it vanishes; hence we should expect the polynomial to be of degree 
less than six. 
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Since we have taken a=0 and A=1, Eqs. (32) and (33) will simplify to 
D(d) = pot piri + port? + + part + por’, 


> A*D(0)b;(s). 


s=t 


If we now write the differences in a column, such that they are spaced the same as 
the elements of Table I, we can readily calculate the coefficients p;. Placing this col- 
umn beside the first column of the table, so that AD(0) is opposite b,(1), and A'D(0) 
is opposite b,(5), we multiply across and add. This gives p;. (If we were dealing with 
a general case in which h#1, we should have to divide by h.) Moving to the next 
column, keeping the same relative vertical position, the operation is repeated to get pz. 
(In the general case we would divide by h?.) Blanks in the table are treated as zeros. 
The calculation runs as follows (the individual products would not normally be writ- 
ten down, but merely accumulated in the calculating machine): 


bi(s) A'D(0) [b:(s)A*D(0)] be(s) A'D(0) 
1 1.00000 — 36 — 36 — 36 0 
2 —0.50000 —476 238 0.50000 —476 — 238 
3 0.33333 —2628 — 876 —0.50000 —2628 1314 
4 -—0.25000  —4968 1242 0.45833 —4968 —2277 
5 0.20000 —2880 — 576 —0.41667 —2880 1200 
Pi =-—1 
Similarly p;= —36, p4=33, ps= —24, and p»=D(0)=20. Hence D(A) =20—8A—X? 


— 36d? + — 24K5. 
c) Method of undetermined coefficients. We assume an expansion of the form 


D(d) = pot prt par? +--+ + pr 


If we assign r+1 values to A, we will obtain r+1 simultaneous linear equations for 
the determination of the p’s: (If we use \ =0 as one value, D(0) = po, so we need solve 
only r equations.) 


This set of equations can then be solved by Aitken’s method [28]. 


™ : In addition to the number of operations (34) for the evaluation of the numerical py 
ee determinants, the method of undetermined coefficients will require 

+ 3) M-D, 

3r(r? — 1) A-S, 


to solve the simultaneous linear equations. This makes a total for this method of 


+ 3) + mr? + (r + 1)(3)(m* + 2m — 3) (36) 


m 
6 


dr(r? — 1) + mr? + (r + 1) (~) (m — 1)(2m — 1) A-S, 


It is apparent that the interpolation formula method requires fewer operations 


= 
| 
| 
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than the method of undetermined coefficients, but if irregularly spaced values of the 
determinants are already at hand, it is better to finish the expansion by the solution 
of the simultaneous equations than to have to evaluate several new numerical deter- 
minants. 

In the special case in which the interval between successive values of X is unity, 
the set of simultaneous equations can be written in the form 


r1 0 0---O07 Pp] 


1 3:9---3 ps D(3) 


Since the matrix A is independent of the numerical values of D(A), we can calculate 
the reciprocal matrix A~! once for all, and then the solution of such a set of equations 
reduces to the formation of the matrix product 


{po Pu Pe} = D(1), , D()}. 

The matr’x multiplicaton on the right will require 

r(r + 1) M-D, 

r? A-S. 
(The first row of A~' will be [1 00 - - - 0], since p»>=D(0).) Comparing these figures 
with the number of operations required for the interpolation formula method, it is 
seen that the latter method will require the same number of additions and subtrac- 
tions, but (4r)(r+1) fewer multiplications and divisions. Consequently it has not 
seemed worth while to tabulate the matrices A~' in this paper, especially when we 
take into consideration the other advantages of the interpolation method, such as the 


information obtainable from the difference table. 
Example. Let us consider thé.determinant 


= 
—2+ rA—A?+ 4% 44+ 2A — 2d? — BA? 
We have 

D(0) = 20 = po, 
D(1) = —16=pot Pit pot Ps + Pat ps, 
D(2) = — 528 = pot 2pit 4p2+ 8ps+ 16f4+ 32ps+ 
D(3) = — 4144 = pot 3pit 22Zpst+ 243p5+ 7729p, 
D(4) = — 18460 = 4pi + 16f2+ 256p4 + 4096p., 
D(5) = — 58920 = pot Spi + 25p2 + 125p3 + 625p4 + 3125p5 + 15625 p¢, 


D(6) = — 151696 = po + 6p1 + 36p2 + 216ps + 1296p, + 7776s + 46656 ps. 


The solution of this set of equations will give the same results as were obtained by 
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use of the interpolation formula, but will require somewhat more labor. ~¢ vanishes, 
as can easily be seen by substituting the column of constants for the ps column in the 
determinant of the coefficients, yielding a vanishing determinant and indicating that 
the polynomial is of fifth degree. 

5. A comparison of methods for obtaining the polynomial expansion. -a) Methods 
applicable to the case |_A _ | =(). In choosing the best method for the expansion of aX 
determinant into polynomial form, it is necessary to consider direct expansion for 
the lower orders. The number of operations required to reach the polynomial form by 
successive expansion in terms of minors of a row or column is given by the recursion 


relations 
M(n) = nM(n — 1) + n(n — 1), M(2) = 2, 


A(n) = nA(n — 1) + 2(n — 1), A(2) = 2, 


(37) 


where M(n) represents the number of multiplications required to expand an mth order 
determinant of the form |.A —Th| , and A(m) the number of additions and subtractions 
for such an expansion. 

We must also consider the application to this case of the interpolation formula 
method (§4b) and the method of undetermined coefficients (§4c). Since the unknown 
appears only on the principal diagonal with unit coefficients, the evaluation of the 
elements of a determinant of mth order will require only n subtractions. Consequently 
the evaluation of the »+1 numerical determinants necessary for either of these meth- 
ods will require (assuming \=0 as one value used) 


(1/3)(n? — 1)(n? + n + 3) M-D, 
(1/6)n(n? — 1)(2n — 1) + n? A-S. 
TABLE II 
Eq. | Method - Mult. and Div. | Add. and Sub. 
| Direct M(n)=nM(n—1)+n(n—1) | 
Expansion M(2) =2 | AQ)=2 

38 | Interpolation Formula 
39 | Undetermined Coefficients (1/6) | 
7 | Leverrier +3) 
11 | Krylov | | (n/6)(n—1)(2n?+7n—1) 
10 | Modified Krylov | G/2)n*(n+1) | (n/2)(n—1)(3n+1) 
12 | Danielewsky | (n*—2)(n—1) | m(n—1)* 
13 | Modified Danielewsky | (n—1)(n*-+n—1) | n(n—1)? 
14 | Reiersgl | (5)(2") | (4)(2") — (n° +-3n+4) 
20 | Samuelson | | 


| 
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For the case in which h=1, the additional operations required to form the difference 
table and the coefficients by the interpolation formula method run the total to 


(1/6)(2n4 + 2n* + 7n? + n — 6) M-D, (38) 
(n/6)(2n* — n? + 10n + 1) A-S. 

The totals for the method of undetermined coefficients are 
(1/6)(2n4 + 5n* + 13n? — 2n — 6) M-D, (39) 
(n/3)(n* + n? + 2n — 1) A-S. 


For ease of comparison, the expressions for the number of operations required for 
each of the methods discussed are placed together in Table II. The relative efficiencies 
of the methods can be seen better from Table III, which gives the actual number of 
operations required by each method for several orders of determinants. 


TABLE III 
Order |} 2 3 4 ae 7 9 
Method | M-D| A-S || M-D| A-S|| M-D! A-S M-D as||mM-p| as || mp | As 
Direct 
Expansion 2} 2 12] 60! 46 || 320 | 238 || 13692 | 10078 || 986400 | 725758 
Interp. | 
Formula | 12) 411 |} 46| 38 |) 125 | 102 |] 279 | 230 || 972| 826) 2525 | 2202 
Undet. Coeffi. | 19 | 10 | 67| 41 || 171 | 116 || 364 | 265 || 1189 | 945 | 2966| 2481 
Leverrier | 6 | 3 41 | 27 || 153 | 114 || 414 | 330 || 1791 | 1533 || 5228] 4644 
| 

Krylov | 17 | 7 || 67 | 38 || 179 | 118 | 389 | 280 || 1287 | 1022 || 3209] 2688 
Modified | | | 
Krylov 18 | 7 || s4| 30 || 120| 78 || 225/160 || sss} 462 1215 | 1008 
Deniclewsky | 2 | 2| 14] 12 || 42] 36 | 92 | so || 252 632 576 
Modified | 
Danielewsky || 5 | 22] 116| 330| 252 712 576 

| 
Reiersgl 3 | 10 43 32 110! 84 | 558 | 438 2438 | 1936 
Samuelson || 4 | 2 || 15 || 48 || 127 | 107 | 369 | 327° 795 723 


Direct expansion proves to be as fast as any method for the second and third 
order ases. Even in the fourth order case none of the methods gives a sufficient 
saving over direct expansion to warrant learning a new technique if only a few equa- 
tions are to be solved. Danielewsky’s method requires the fewest operations from the 
fifth order up, but it is really harder to use than any of the last three methods, which 
are about equal in the fifth order case. Above the fifth order, the two most efficient 
methods are Samuelson’s and the Modified Danielewsky. The one to be used will de- 
pend a great deal on the habits of the computer. 
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pre If we have already started solving the secular equation by the matrix iteration 

4 ‘ method, we will have computed all or part of the sequence C(k) of Krylov’s method, 

2 and it will usually be quicker to complete the polynomial expansion by that method. } 
é . If we have already evaluated D(A) for several values of \, as we might do in hunt- 

4 : ing for a root by the method of false position, it would be preferable to finish by the 

{ interpolation formula method or the method of undetermined coefficients, depending ' 
me. on whether the successive values of \ are uniformly spaced or not. 

; . If a machine were available on which matrices could be multiplied with ease, 


Leverrier’s method would be useful, since the set of simultaneous equations for the 
coefficients is so simple. 


al b) Methods applicable to the case |.A — By| =0. The number of operations required 
» se for direct expansion of the determinant | A — By| is given by the recursion relations 
o M(n) = nM(n — 1) + 2nt M(1) = 0, 

40 
_ A(n) = nA(n — 1) + (m — 1)(2n + 1) A(1) = 0. wi 


The various methods applicable’ to this case are compared in Tables IV and V. 


TABLE 1V 
3 Eq. Method | Mult. and Div. | Add. and Sub. 
: 40 | Direct Expansion M(n) =nM(n—1)+2n? | A(n) =nA(n—1)+(n—1)(2n+1) 
35 Interpolation Formula | (1/6) (2n4+8n?+7n?+n—6) | 


36 | Undetermined Coefficients | (n/3)(n?+-4n? —n —1) 


(n/2)(n—1)(Sn—3) 


22 | Reciprocation | 


23 | | | 


TABLE V 


Order | 2 | 3 a | s | 7 9 
| | | | 
é | | | } 
. Method | M-D) A-S| M-D| A-S|| M-D| A-S || M-D] A-S || M-D | AS |} M-D | A-S 
Direct | | | f 
Expansion || 8 | 5 || 42 | 29 || 200 | 143 |/1050 | 759 | 44702 | 32423 || 3219858 | 2335679 
| 
Formula | 20 | 15 || 73 | 56 || 189 | 150 | 404 | 330 | 1315 | 1120 3254 | 2850 
Undet. | | | | 
* Coeff. | 27 | 14 | 94 | 59 || 235 | 164 || 489 | 365 || 1532 | 1239 3695 3129 
Reciprocation || 15 | 7 | 58 | 36 | 145 | 102 | 291 | 220 || 672 1765 | 1512 
Danielewsky- | | 
Masuyama || 10| 6 || 42 | 30 || 115 | 89 255 | 207 || 875| 751 2250 1994 
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Again direct expansion is to be preferred for second and third order cases. The 
Danielewsky-Masuyama method requires the fewest operations in the fourth and 
fifth orders, but because of the large amount of writing necessary with this method, 
the author prefers the method of reciprocation for all cases from the fourth order up. 

In case the determinant | B| vanishes, the interpolation formula method is to be 
preferred except in the case in which one has already obtained several values of 
| A —BX| for unequally spaced values of d, in which case the method of undetermined 
coefficients is preferable. 

c) Methods applicable to the case |A +Br+0r?3| . If m is the order of the determi- 
nant, direct expansion will require the number of operations given by the recursion 
relations 

M(m) = mM(m — 1) + 3m(2m — 1), M(1i) = 0, 
A(m) = mA(m — 1) + (m — 1)(6m + 1), A(1) = 0. 


The various methods applicable to this case are compared in Tables VI and VII. 


(41) 


TABLE VI 
Eq. Method Mult. and Div. Add. and Sub. 
41 Direct M(m) =mM(m—1)+3m(2m—1) A(m) =mA(m—1)+(m—1)(6m+1) 
Expansion M(1)=0 A(1)=0 
30 Transformation 9m* —2m? —6m+2 9m* —13m?+6m+2 
35 Interpolation 
Formula (1/3) (2m*+-13m?+ 10m? —m —3) (m/6) (4m*+20m?+23m +1) 
36 Undetermined 
Coefficients (1/3) (2m*+-25m*+22m*? —4m —3) (m/6) (4m +44m?—m 
TABLE VII 
Order | .' Ess 4 5 6 7 
Method | M-D | AS | M-D) A-S || M-p| A-S|| M-D| A-s || | || M-D | as 
Direct 
Expansion ] 18 13 99 77 || 480 | 383 || 2535 | 2039 || 15408 | 12419 || 108129 | 87191 
Trans. | 54 34 || 209 | 146 || 522 | 394 || 1047 832 1834 1514 2949 2494 
Interp. | 
Formula | 57 | 53 || 199 | 179 || 499 | 446 || 1039 | 930 1917 1723 3247 | 2933 
Undet 
Coeff. 103 67 || 340 | 248 || 815 | 634 || 1634 | 1325 2919 2437 4808 4102 


In this case, direct expansion requires the smallest number of operations through 
the fourth order, the interpolation method the smallest number for the fifth order, 
and the method of transformation to the diagonal form the smallest number for all 
higher orders. The interpolation formula method is convenient, however, even when 
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it requires some extra operations, as it makes possible a preliminary plot of the func- 
tion. The interpolation formula method is the best method above the fourth order 
when the determinant | C| vanishes, as the difference table gives the order of poly- 
nomial to be expected. 

The method of undetermined coefficients is useful primarily when several values 
of the determinant have already been calculated for unevenly spaced values of X. 

6. Errors. Errors of considerable magnitude can readily occur in the various proc- 
esses described in this paper, due both to errors in the original data and errors caused 
by rounding off numbers. The study of errors arising in the evaluation of determi- 
nants, solution of simultaneous linear equations, iterative methods, etc., is far from 
complete. A discussion of these errors is beyond the scope of this paper, but the inter- 
ested computer should consult references [29-34]. 


REFERENCES 


1. JAMEs and CootipcE, J. Chem. Phys. 1, 825 (1933). 
2. DuNncAN and Cotxar, Phil. Mag. 17, 865 (1934). 
3. FRAZER, DUNCAN and CoLar, Elementary matrices, Macmillan, 1938, p. 133. 
4, A. C. AITKEN, Proc. Roy. Soc. Edinburgh, 57, 269 (1937). 
5. T. C. Fry, private communication. 
6. M. Masuyama, Proc. Phys. Math. Soc. Japan (3), 21, 100 (1939). 
7. A. N. Krytov, Bull. de l’Acad. Sci. U.R.S.S. (4), 7, 491 (1931). 
8. A. DANIELEWSKY, Trans. Mosc. Math. Soc. 44, 169 (1937). 
9. FrAzER, DUNCAN and CoLtar, loc. cit., 141-142. 
10. B. Hicks, Thesis, Calif. Inst. of Tech. (1939). 
11. B. Hicks, J. Chem. Phys. 8, 569 (1940). 
12. P. Horst, Ann. Math. Stat. 6, 83 (1935). 
13. P. A. SAMUELSON, Ann. Math. Stat. 13, 424 (1942). 
14. O. Rerers@x, Ann. Math. Stat. 11, 193 (1941). 
15. F. Cu1d, Mémoire sur les fonctions connues sous le nom de résultants ou de déterminants, Turin, 
1853. Or, E. T. Whittaker and G. Robinson, The calculus of observations, 2nd. ed., Blackie, 1929, p. 71. 
16. A. C. AITKEN, Proc. Roy. Soc. Edinburgh 46, 289 (1926). 
17. A. C. AlTKEN, Proc. Roy. Soc. Edinburgh 51, 80 (1931). 
18. L. A. Pires, Jour. Franklin Inst. 225, 437 (1938). 
19. U.J. J. LEVerrtErR, J. de Math. s.1, 5, 230 (1840). 
20. Dickson, First course in the theory of equations, Wiley, 1922, p. 134. 
21. Mutr and Metz_er, Theory of determinants, Albany, N. Y., 1930, p. 606. 
22. FRAZER, DUNCAN and COLLar, loc. cit., p. 141. 
23. BOcHER, Introduction to higher algebra, Macmillan, 1927, p. 296. 
24, FRAZER, DUNCAN and COLLar, loc. cit., p. 90. 
25. Ibid., p. 162. 
26. WHITTAKER and RoBINSON, loc. cit., p. 10. 
27. Ibid., p. 28. 
28. A. C. AITKEN, Proc. Roy. Soc. Edinburgh 57, 172 (1937). 
29. I. M. H. EtHERINGTON, Proc. Edinburgh Math. Soc. (2) 3, 107 (1932). 
30. F. R. Moutton, Amer. Math. Monthly 20, 242 (1913). 
31. L. B. TucKERMAN, Annals Math. Stat. 12, 307 (1941). 
32. P. G. Hoe, Ann. Math. Stat. 11, 58 (1940). 
33. A. T. Lonsetu, Ann. Math. Stat. 13, 332 (1942). 
34. Hore Ann. Math. Stat. 14, 1 (1943). 


: 

if 

| 

4 

4 

— 9 

aa ly 

Ps. 


THE PROBLEM OF SAINT VENANT FOR A 
CYLINDER WITH FREE SIDES* 


BY 


J. L. SYNGE 
The Ohio State University 


1. Introduction. Consider a cylinder composed of homogeneous isotropic elastic 
material (Fig. 1). The cross section is arbitrary; it may be simply or multiply con- 
nected. Body force (such as gravity) is assumed to be absent, and the sides are free. 
To the ends we apply any loadings 
which satisfy the conditions of x 
statical equilibrium. As a result ' 
the cylinder undergoes a small de- 
formation. We ask: what is the 
stress and what is the displace- 
ment throughout the cylinder? 

Here we have a well-formu- 
lated problem in the theory of O 
elasticity, and one of the most im- \/ 
portant from a practical stand- 
point. It may be called the problem 
of Saint Venant. It includes as spe- 
cial cases the problems of tension, 
bending by couples, torsion and 
flexure. 

Saint Venant' ingeniously side-stepped the major difficulty of the problem by 
substituting a simplified problem. Of the six components of stress (the unknowns of 
the original problem) he set three equal to zero, and relaxed the boundary conditions 
on the ends of the cylinder. In his solution the terminal loadings are not arbitrarily 
distributed; only the total load can be arbitrarily assigned. In his solution of the ten- 
sion problem, for example, the load must be uniformly distributed over the cross 
section. 

In tension or bending by couples, Saint Venant’s solution is mathematically 
trivial. For torsion or flexure, the mathematical high spot is a problem of Dirichlet 
or Neumann (determination of a function harmonic in the cross section, the values 
of the function, or of its normal derivative, being assigned on the boundary). Mathe- 
maticians have busied themselves through the years with general and particular 
methods of solution. The torsion-flexure problem has been regarded as one of the 
central problems in the theory of elasticity. 

Mathematicians should, however, be reminded that when they speak of the “tor- 
sion-flexure problem,” they have in mind Saint Venant’s type of solution, in which 
the terminal loads are not arbitrarily distributed. Relatively little attention has been 
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1 B. de Saint Venant, Journal de Math., 1, 89-189 (1856). 
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paid to more general types of solution. It is good, I think, to draw attention to the 
work that has been done, and to indicate possible methods of attack. 

Section 2 gives the general mathematical statement of the problem of Saint 
Venant, and of the relaxed problems associated with it. Section 3 contains an old 
friend in new clothes—the Saint Venant solution in tensor notation. I believe that 
is worth while to give such a presentation, avoiding as far as possible the customary 
special choices of axes and particular integrals, in order to reveal the true mathe- 
matical structure of the solution. In Section 4 we pass beyond the Saint Venant solu- 
tion, and set up the basic eigen-value problem associated with the exponential type 
of solution. In Section 5 we link this formulation with the solution given by Dougall 
for a cylinder of circular section. Section 6 contains some questions. 

2. The problem of Saint Venant and the relaxed problems. Latin suffixes will have 
the range 0, 1, 2, and Greek suffixes the range 1, 2, with the usual summation conven- 
tion for repeated suffixes. Let x; be rectangular Cartesian coordinates, with the axis 
of x» in the direction of the generators of the cylinder. The position of the origin and 
the directions of the axes of x; and xz remain arbitrary. Let u; be the displacement and 
Ei; (=£,) the reduced stress, i.e., the stress divided by Young’s modulus. We have 
the basic stress-strain relations 


+ = (1 + o) Bij — (2.1) 


or equivalently 


2(1 + o)(1 20) E;; = 205; + (1 20) + Uj, (2.2) 


Here a is Poisson’s ratio (a constant which may, in theory, take any value in the range 
—1<o<4), and the comma denotes partial differentiation (f,;=0f/0x;); 5;; is the 
Kronecker delta. 

In any problem in elasticity, we have a choice between two methods: we can work 
with the displacement uw; or with the stress £;;. When the boundary conditions are 
given in terms of stress (as they are in the problem of Saint Venant), the relative ad- 
vantages may be set down as follows: 

Displacement method: Simple p.d.e. and complicated b.c. 

Stress method: Complicated p.d.e. and simple b.c. 

Here, and throughout, “p.d.e.” means “partial differential equations,” and “b.c.” 
means “boundary conditions.” 

The two rival formulations of the problem of Saint Venant are set out below in 
(2.3) and (2.4). In each case, (a) contains the p.d.e., (b) contains the b.c. on the free 
sides (m, are the direction cosines of the normal), and (c) contains the b.c. on the 
ends. 7», T are the components of the assigned stress. The symbol A; is the 3-dimen- 
sional Laplacian differential operator (A3=0?/0x;0x;). 


Saint Venant problem in terms of displacement: 
(1 20) Agu; + 0; (2.3a) 
(ug,0 + = 0, + (1 — 20)(Us,0 + = 0; (2.3b) 
OUuUk,k + (1 20) uo,0 ="(1 2¢)To, U0,a + 2(1 + o)T a. (2. 3c) 
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Saint Venant problem in terms of stress: 


E;;,; = 9, (1 + o)AsEi; + = 0, (2.4a) 
Engng = 0, E.gns = 0, (2.4b) 
Eo = To, Eoa = Ta. (2.4c) 


We note that in (2.3) there are 3 unknowns, 3 p.d.e., and 3 b.c. In (2.4) there are 6 un- 
knowns, 9 p.d.e., and 3 b.c.; but between the 9 p.d.e. there exist 3 differential identi- 
ties. 

Let us cross out the (c) equations in (2.3) and (2.4), that is, drop the conditions 
on the ends of the cylinder. We have then alternative statements of what may be 
called the relaxed problem of Saint Venant. The unknowns are of course underdeter- 
mined in the relaxed problem, but the system is now linear and homogeneous, so that 
solutions may be superimposed. We may hope that, by superimposing solutions of the 
relaxed problem, we may succeed in satisfying the b.c. (c), either accurately or ap- 
proximately. 

In spite of the underdetermination, the relaxed problem is too complicated to yield 
to mere guessing, except for one very simple solution: Zoo =const., Eo. =0, Eag=0. 
Indeed, it is usually more difficult to deal with an indeterminate problem than with a 
determinate one, and so we impose auxiliary conditions to replace the b.c. (c). 

The best-known auxiliary condition is the so-called hypothesis of Saint Venant: 


Ea = 0. (2.5) 


This leads to the Saint Venant solution, which will be discussed in Sect. 3. Another 
auxiliary condition, which we may call the exponential condition, is 


Ey; = x2), (k #0). (2.6) 
This is suggested by the fact that if Z;;=f,; satisfy the relaxed problem, then so also 
do E;;=fi;,0. The problem of determining solutions subject to (2.6) will be discussed 


in Sect. 4. 
Dougall? has used the expression permanent free modes for solutions of the relaxed 


problem with (2.5) imposed, and transitory free modes for solutions of the relaxed prob- 
lem with (2.6) imposed. (This terminology is suggested by the theory of vibrations, 
but of course our problem is statical.) According to Dougall, these two types of solu- 
tions are fundamental, in the sense that any solution of the relaxed problem is a 
linear combination of them. 

3. Saint Venant’s solution. Any treatise on elasticity contains a treatment of 
Saint Venant’s solution of the relaxed problem, usually broken up into the problems 
of tension, bending; torsion, and flexure* for pedagogic reasons. Clebsch® appears to 
have been the first to see that Saint Venant’s solution follows logically from (2.5). 
Marcolongo® has given an elegant general treatment, using displacement as funda- 
mental. In the following treatment, stress is used and the axes 0x;x2 remain completely 
general. 


2 J. Dougall, Trans. Roy. Soc. Edinburgh, 49, 895-978 (1913). See also Proc. Fifth International 
Congress of Mathematicians, 2 (Cambridge 1913), 328-340. 

3A. E. H. Love, Mathematical theory of elasticity (Cambridge, 1934), pp. 329-334. 

41. S. Sokolnikoff, Elasticity, Bréwn University, 1941, pp. 193-202. 

5 A. Clebsch, Theorie der Elasticitat fester Kérper, Leipzig, 1862. 

6 R. Marcolongo, Teoria matematica dello equilibrio dei corpi elastici, Milano, 1904, pp. 296-310. 
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We substitute the auxiliary conditions (2.5)-in (2.4a, b), and obtain 
Eoo,o + Eog.s = 0, (1 + o)AsEoo + Eoo,oo = 0, 
E.0,0 = 0, (1 + o)AsEao + Eoo,a0 = 9, (3. 1a) 


Eoo,ap = 0; 
Eogng = 0. (3. 1b) 


It follows at once that E,o are independent of x, and that 
Eoo = Xo(A A) Bsxs B, (3.2) 
where Az, A, Bs, B are six constants, at present arbitrary. The remaining equations 
in (3.1a) are equivalent to 
AEa = (1 + Es0,8 = Agxs + (3.3) 
where A is the plane Laplacian (A = 02/0xg0xg). These three equations are to be solved 
for the two unknowns E,o, with the b.c. (3.1b). The problem is a plane problem, the 


domain being the section of the cylinder. 
Given the vector E.o, there exist invariants ¢, Y, such that 


Ea = $,a — €ayW,7 (€11 = €22 = 0, €12 = — = 1). (3.4) 
Substitution in (3.3) leads at once to 
Ad = + A, (1+ = +C, (3. 5a) 
where C is another arbitrary constant. (Note that €ag€a7 = 5sy.) The b.c. is 
— = 0 or d¢/dn — = 0, (3.5b) 


where 9m is a element of the normal m,, drawn out of the material and to the right of 

the element 0s of the bounding curve 

C (Fig. 2). On integrating (3.5b) around 

the bounding curve, and using (3.5a), 
we find that 


A = — (3.6) 


where # are the coordinates of the 
centroid of the section. Thus there 
are only six arbitrary constants, Ag, 
Bz, B, C. 

Remembering the Riemann-Cau- 
chy relations, we note that ¢, W are 
indeterminate to the extent of adding 
to 6+ an arbitrary analytic func- 
tion of x; 

Let be any single-valued 


Fic. 2, particular solutions of (3.5a). Let us 
define by 
(3.7) 


so that @, YW are harmonic. Let 2 be the harmonic conjugate of V, so that 
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Qe = = ay. (3.8) 
If P is defined by 
P=4-Q, (3.9) 
then (3.4) reads 
and the b.c. (3.5b) reads 
OP/dn = — 06 /dn + dy /ds. (3.11) 


Thus the relaxed problem of Saint Venant reduces to a problem of Neumann—to 
find a harmonic function P under the b.c. (3.11). The solution must be single-valued 
in order that the displacement may be single-valued. 

This presentation of the problem has a fictitious simplicity; the constants As, C 
are wrapped up in the Neumann problem. We may split up the problem into problems 
involving only the geometry of the section and Poisson’s ratio. To do this, we write 


P = AgPs + (1 + 


(3.12) 
Here ¢$”, vf”, x are single-valued particular solutions of 
= — (1+ = -Ax” = 1, (3. 13) 


and Ps, Q are single-valued harmonic functions satisfying the b.c. 
= — + 9Q/dn = dx /ds. (3. 14) 


The determination of the vector Ps constitutes the flexure problem, and the determi- 
nation of Q the torsion problem. 
Since =8x., A(r?)=4, where r?=xgxg, particular solutions of (3.13) are 
given by 


(1 + = orreg,x,/8, x) = 


The determination of the displacement in terms of (3.2) and (3.10) in this general 
notation is interesting, because it reveals why the components u, are independent of 
the solution of the Neumann problem (3.11). Using (2.5) in (2.1), we have 


= Evo, Ua,0 + = 2(1 + Ea, 


(3.15) 


(3.16) 
UB + 206 agE oo. 
By (3.2), the first two equations give at once 
= — A B B 
uo $x9(Agxg + A) + x0(Baxs + B) + fo (3.17 


Ua = 4x3A. = Xofo,a + + a) ao + fa 


where fo, f. are unknown functions of x;, x2. When we substitute this expression for u. 
in the last of (3.16), we get 

fo.as = (1 + + Epo,a) — + A), (3.18) 

Sea t+ fap = — 206 ap(Byxy + B). (3.19) 


a 
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Now, by (3.10), 

E206 — Epo,a = — + ep, = — eaphy™, (3.20) 
and so, since y“ satisfies (3.5a), 


(1 + Es0,a) €ap (CE + (3.21) 


We now use this equation to sbustitute for Zgo,4 in (3.18). When we do so, the linear 
expression on the right must be a partial derivative with respect to xs. Now if Fas... 
is a single-valued homogeneous function of degree m in x;, x2, it follows from Euler’s 


theorem that 


1 
= — Fag... + (3.22) 
n 


where K is a constant. Hence we deduce from (3.18) 


= 21+ a) Eao + + C) TXa(ZA yXy + A) + he, (3.23) 


where h, is constant. Now substitute for EZ.» from (3.10) and integrate; this gives 


fo = 21 + o)P + + + [[— 201 + + yt, + C) 
— + A) |dxa + hata 


(3.24) 


where / is constant. By virtue of the p.d.e. (3.5a) satisfied by y, this integral has 
the same value for reconcilable paths. But we have no guarantee that it has the same 
value for irreconcilable paths in the case of a multiply connected section. To secure 
a single valued displacement, we should choose for y“ a function which satisfies the 
p.d.e. (3.5a) throughout the whole interior of the outer boundary. We may, for ex- 
ample, use the expression given in (3.12), (3.15). 

The equations (3.19) determine f. to within a plane rigid body displacement, 


and so 
fa = Bar? — — + heayXy + ka, (3.25) 


where &,, & are constants. 

To find the displacement, we substitute for fo from (3.24) in the first of (3.17), 
and so obtain up. It involves the solution of the Neumann problem, since P appears 
in (3.24). To find u,, we substitute in the second of (3.17). There are two substitutions, 
one for fo,2—2(1+0¢)Eao from (3.23), and the other for f, from (3.25). We see that P 
does not occur, and thus we verify the well known fact that the lateral displacement 
is independent of the solution of the Neumann problem. 

4, The exponential type of solution. Let us now investigate the solution of the 
relaxed problem with the auxiliary condition (2.6). Since stress determines displace- 
ment to within a rigid body displacement (which we shall omit), we may write for the 


corresponding displacement 


Once more we have a choice between two methods—displacement and stress— 


and we shall choose displacement. 
When we substitute from (4.1) in (2.3a, b), the equations reduce to 
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(4 + k?)V = 0, (4 + k*) Vie (4. 2a) 
20V na a (06,0 + Ng = 0, (i 2c) V + k°vgng — = 0. (4. 2b) 


Here V is an auxiliary variable, given by 
(1 — = hoo + 0,8; (4.3) 


it is the dilatation, to within a factor. 

Inspection of (4.2) shows that we have before us an eigen-value problem of con- 
siderable complexity. The system will (presumably) be consistent only for certain 
values of k. There is no objection to complex eigen-values, with complex solutions 
for V, v«, and the corresponding stress. Denoting complex conjugates by bars, we 
should take in such cases for the real displacement and stress 


+ + e*°F;;). (4.4) 


If & is an eigen-value of (4.2), so also are —k and +&. In fact, the eigen-values occur 
in sets of two if they are real or purely imaginary, and in sets of four if complex. 

By a simple and ingenious argument, Dougall? has shown that the system (4.2) 
has no purely imaginary eigen-values. A purely imaginary k implies a periodic dis- 
tribution of displacement and stress. Consider the energy in a length of cylinder equal 
to this period. It is equal to the work done by the terminal stress in passing from the 
natural state to the strained state. But, from the periodicity, this is zero. Hence, the 
energy of a strained state is zero, which is contrary to a basic postulate of elasticity. 
Hence there can be no purely imaginary eigen-value k. It should be added that we can- 
not assert this if o is arbitrary. It is necessarily true only if strain-energy is positive- 
definite, i.e., if —1<o0<}. 

In many problems in applied mathematics, harmonic functions play a funda- 
mental role. In our problem (4.2) that role is taken over by plane wave functions, 
where by a wave function f(x:, x2) we mean a solution of (A+£?)f=0. We note that V 
is a wave function, but v, is not. 

We can, however, easily reduce the unknowns to wave functions by writing 


Vq = Wa + We, (4.5) 
where w* is any particular solution of 
(A + k*)wt = — (4.6) 


Let us not tie ourselves down to any definite particular solution. Two, however, seem 


to be particularly interesting: 
we = — (4.7) 


we = — Wa (A+ k*)W = V. (4.8) 

Our problem (4.2) may now be stated as follows: 
(A + k)V = 0, (A + = 0, (4. 9a) 
ina + (Wha + We + + = 0, 


(1 — 20)V png + — + k*weng — = 0. 


(4. 9b) 


In simplifying the p.d.e., we have complicated the b.c. 
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In dealing with Saint Venant’s solution, we found it advisable to pass in (3.4) 
from a vector to two invariants. The same procedure will be pursued here. We shall 
make use of the following fact: Given a vector wave function Wa, there exist wave functions 


such that 
Wa = — (4.10) 


and >, are unique. 
To prove this, we have merely to put 


It is easy to verify that these satisfy (4.10), by virtue of the wave character of w,. As 
for uniqueness, we can add to ¢, ¥ only functions ¢®, y™ such that ¢+i~™ is 
an analytic function of x1+7x2. Thus ¢®, y™ are simultaneously wave functions and 
harmonic functions; therefore they vanish. 

Let us now transform the b.c. (4.9b). Let C (Fig. 2) be the boundary curve of the 
cross section. Let ¢, =dx./ds. Then 


Ne = €apls, 


(4.12) 
dn,/ds = kta, dt,/ds = — 


where x is the curvature of C, positive if C bends away from mg. 
We may replace the first of (4.9b) by two invariant equations, obtained by multi- 
plying by ma, ta, respectively. Thus 


oV + + = 0, 


(4.13) 
+> + (Wp,a + Wapltang = 0. 
From (4.10) we have 
+ = 2, apnang — aptans, (4.14) 
+ = 26, astang + — WV, aptals. 
The following identities are easily verified: 
0°/dn? = ~(*) + Kb, at 
as \an (4.15) 
Ad = $,asnang + ¢$,aptals, 0°/ds? = apslalg — KP,aNa- 
Hence 
(4.16) 
( + )t 2(< + + 


Thus our eigen-value problem may be stated in the following final general form: 
Given a domain in the (x1, x2) plane, bounded by a curve C, we seek eigen-values k and 
eigen-functions V, , W to satisfy the p.d.e. 


(A+ k)V = 0, (A + k*)¢ = 0, (A + k*)y = 0, (4.17a) 
and the b.c. on C 
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2 
oV + - + k? + -(2 = 0, 
ds? on Os On Os 


(was + + 4k 0 ( ) 
(1 + 20) — + — + 2k? Boe 
on on Os 


Here w,* is any particular solution of (4.6). 
If we choose w,* as in (4.7), we may substitute in (4.17b): 


‘ 1 1 av 
oV + Wa,pnang= —— (1-— — 
2 4 on 


— —— — ——— —(r) 


on On 4 on 2 On 


On the other hand, if we choose w.* as in (4.8), we may substitute for these three 
quantities, respectively, 
ow 
oV —a°W/ant, — (— =) ¥, —- 6.9 
Os On Os On on 
If we choose w,* as in (4.7), we have the following expressions for the displacement 
in terms of V,¢, ¥, k 
Uo = Ua = e**,, 
= 2(1 o)V + $z.V.. + (4.20) 
It takes only a moment to verify that these expressions satisfy (4.2a). These formulas 
(in a slightly more general form) are the key to Dougall’s treatment of the circular 
cylinder.? However, he gives no indication as to how to he obtained them, or whether 


they are a general representation of the exponential type of displacement. I have 
shown that they are in fact a general representation. Given v; and k, then V, ¢, ¥ 


are uniquely determined. 
The stress corresponding to (4.20) is 


E,;=e**F;;, Fao= — Foo= k-*F 28,08; 
2(1+¢)Fas= p+ xpV a) +26, (Cav. 7a), (4.21) 
2(1+6)Foo=2(2—0)V a 


5. The case of the circular section. The circular section has been dealt with so 
thoroughly by Dougall? that there is little more to be said about it. However, it is 
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In dealing with Saint Venant’s solution, we found it advisable to pass in (3.4) 
from a vector to two invariants. The same procedure will be pursued here. We shall 
make use of the following fact: Given a vector wave function Wa, there exist wave functions 


such that 
Wa = 9,a — (4. 10) 


and are unique. 
To prove this, we have merely to put 


= — k-*e,,w,,». (4.11) 


It is easy to verify that these satisfy (4.10), by virtue of the wave character of w,. As 
for uniqueness, we can add to ¢, only functions such that is 
an analytic function of x,;+7x.. Thus ¢”, y™ are simultaneously wave functions and 


harmonic functions; therefore they vanish. 
Let us now transform the b.c. (4.9b). Let C (Fig. 2) be the boundary curve of the 


cross section. Let tg =dx./ds. Then 
Na = €apls, = 
dn,/ds = kta, dt,/ds = — kNg, 


(4.12) 


where x is the curvature of C, positive if C bends away from mq. 
We may replace the first of (4.9b) by two invariant equations, obtained by multi- 
plying by a, ta, respectively. Thus 


oV + we gnang + = 0, 


(4.13) 
(wea + + + Wa,ptans = 0. 
From (4.10) we have 
(ws, a + Wa,p)NaNg (4 14) 
(Ws,a + Wapltats = 26, aptams + — W,aptals. 
The following identities are easily verified: 
0°6/dn? = + Kb, at 
as an aplallg (4.15) 
Ad = apNang + ¢,aptats, = — KP,aNa- 
Hence 
(4.16) 


Thus our eigen-value problem may be stated in the following final general form: 
Given a domain in the (x1, x2) plane, bounded by a curve C, we seek eigen-values k and 
eigen-functions V, , to satisfy the p.d.e. 


(A + k)V = 0, (A + k*)o = 0, (A + k*)y = 0, (4.17a) 
and the b.c. on C 
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on 


00a 3? 
+ Ws, + (= (= + $k? + = 0, (4.17b) 
Os On Os on 
oy 


(1 2c) k? we ng — We pyny + 2k? — — k? — = 0. 
On on Os 


Here w,* is any particular solution of (4.6). 
If we choose w,* as in (4.7), we may substitute in (4.17b): 


2 4 On on 


On On 4 on 2 On 


On the other hand, if we choose w,.* as in (4.8), we may substitute for these three 
quantities, respectively, 
0a OV 
oV —a*W/an?, — (= =) —- 49 
Os On Os on on 
If we choose w,* as in (4.7), we have the following expressions for the displacement 
in terms of V, ¢, y, k: 
Up = e**%y9, Ua = 
kup = 2(1 a)V + $<... + (4.20) 
It takes only a moment to verify that these expressions satisfy (4.2a). These formulas 
(in a slightly more general form) are the key to Dougall’s treatment of the circular 
cylinder.? However, he gives no indication as to how to he obtained them, or whether 


they are a general representation of the exponential type of displacement. I have 
shown that they are in fact a general representation. Given v; and k, then V, ¢, ¥ 


are uniquely determined. 
The stress corresponding to (4.20) is 


E,;=e'F;;, — kF ag, Foo= k-*F 
2(1+¢)Fas= —(1—2¢)bapV —3(xaV p+ 2) +26,08— (Ca, (4.21) 


5. The case of the circular section. The circular section has been dealt with so 
thoroughly by Dougall*® that there is little more to be said about it. However, it is 
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interesting to see how his method of solution connects up with the preceding ap- 
proach. 
Let the boundary be r=a, and let w be the polar angle. Then the following are 
wave functions 
(V, ¥) = (A, B, C)J (5.1) 
where A, B, C are any constants. Each of these functions satisfies on r=a the equa- 


tions 
of Kf af imf of imKf 


= 


on a Os a Os On a? (5.2) 
K = = ka. 


They also satisfy for r=a, 
— (xfs) = (m? — &)f/a. (S.3) 
On 

Substitution in (4.18) and (4.17b) gives 


LiVa? + Lid + Liew = 0, 
= 3(1 — 2c) + 3K, In = —m’+?+K, Loe = im(K — 1), 


(5.4) 
Lo = — im/4, Lu=im(K-1), Ly 
Loy = (1 — o)K + jm? — In = Ké, Loe = — 
The characteristic equation is 
| Liz] = 0, (5.5) 


which is essentially the determinantal equation of Dougall for the determination of &, 


and hence &. 
For modes independent of w, we put m =0. Then (5.5) reduces to 


(& + 2K) [&(K? + #) — 2K*(1 — =0, (5.6) 
and we get the two characteristic equations 
J2(&) = 0, + Jé?) = 2(1 — (5.7) 


The first equation has an infinite sequence of real roots; the second has an infinite 
sequence of complex roots. 

Dougall states (p. 902) that for every m, (5.5) has an infinite number of real roots 
and an infinite number of complex roots (but no purely imaginary roots, as pointed 
out earlier). 

For other work bearing on the circular section, reference may be made to Poch- 
hammer,’ Thomae,® Schiff,* Chree,’° Filon,!® Purser,!7!8 Timpe,!® 


7L. Pochhammer, Journal f. Math., 81, 33-61 (1876). 

8 J. Thomae, Berichte Verh. K. Sachs. Ges. d. Wiss. zu Leipzig, Math. Phys. Cl., 37, 399-418 (1885). 
® Schiff, Journal de Math., 9, 407-424 (1883). 

10 C, Chree, Cambridge Philosophical Transactions, 14, 250-369 (1887). 

11 QO, Tedone, Atti R. Acc. Lincei Rend. Cl. sci. fis. mat. nat., 10, 131-137 (1901). 

12. Tedone, Atti R. Acc. Lincei Rend. Cl. sci. fis. mat. nat., 13;, 232-240 (1904). 

13 OQ. Tedone, Atti R. Acc. Lincei Rend. Cl. sci. fis. mat. nat., 202, 617-622 (1911). 
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Wolf,?° Barton.*! Barton makes use of the general Papcovitch?*-Neuber**.* solution 


of (2.3a), viz., 


1 
= Yi = (xi + 9), (5.8) 
where A;W¥; =0, A3¢ =0. Any solution of (2.3a) may be so expressed, but it seems that 
for present purposes Dougall’s expressions (4.20) are more convenient. 
6. Conclusion. The method of Dougall can be extended to a pipe, i.e., a section 
bounded by two concentric circles. Here are some questions: 
1) Are there any other sections which can be solved by simple extensions of the 
method used for the circle? 
2) Do eigen-values exist under reasonably general assumptions regarding the 
smoothness of the boundary curve? 
3) Is there always a set of real eigen-values, or is that a peculiarity of the circle? 
4) Write k=p+ig, and let m be the least value of | p| in the sequence of eigen- 
values for a given section. Then m?S, where S is the area of the section, depends 


only on the shape of the section. For arbitrary sections, m*S forms a positive se- | 


quence. Is it bounded below, and, if so, what is the lower bound? 

This last question is very interesting from a practical point of view, because | p| 
represents the rate at which end effects decay as we pass along the cylinder. The 
greater | p| , the more rapid the decay. Engineers are worried by end effects, because 
the Saint Venant solution gives no information about them. The assignment of such 
a lower bound might be more valuable than the description of a complicated process 
for the evaluation of eigen-values. 


4 OQ. Tedone, Atti R. Acc. Lincei Rend. Cl. sci. fis. mat. nat., 21,, 384-389 (1912). 
16 QO. Tedone, Encyk. d. math. Wiss., 44, p. 150. 

1 L. N. G. Filon, Phil. Trans. Roy. Soc. London, A 198, 147—233 (1902). 

17 F, Purser, Trans. Roy. Irish Academy, 32 A, 31-60 (1902). 

18 F, Purser, Proc. Roy. Irish Academy, 26 A, 54-60 (1906). 

19 A. Timpe, Mathematische Annalen, 71, 480-509 (1912). 

20K. Wolf, K. Akad. d. Wiss. Wien, Math.-naturwiss. KI., Abt. IIa, 125, 1149-1166 (1916). 
21 M. V. Barton, J. App. Mechanics, 8, A-97-A-104 (1941). 

22 P. F. Papcovitch, Comptes rendus Acad. des Sci. Paris, 195, 513-515 (1932). 

23 H. Neuber, Zeit. angew. Math. u. Mech., 14, 203-212 (1934). 

24H. Neuber, Kerbspannungslehre, Berlin, 1937, 
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POWER SERIES EXPANSIONS OF THE VELOCITY 
POTENTIAL IN COMPRESSIBLE FLOW* 


BY 


G. K. BATCHELOR 
Division of Aeronautics, Council for Scientific and Industrial Research, Australia 


1, The equation for compressible flow. The equation describing the two-dimen- 
sional irrotational flow of a compressible fluid is 


99 i th og? hed lt 1)(q? (1.1) 
or? , or Or Or rd0 
where r, 6 are polar coordinates, M is the Mach number of the flow at some distance 
from the origin where the velocity is assumed to be uniform, y is the ratio of specific 
heats of the fluid and ¢, g are the velocity potential and velocity, respectively, at the 
point 7, 6. All distances are referred to some typical dimension in the field, and veloci- 
ties are referred to the speed of the uniform stream. The boundary conditions to be 
imposed on ¢ are (1) dg/dr—cos 6 as r— ~, where @ is measured from the direction of 
the uniform stream, and (2) the normal derivative of g equals zero at the boundary 

of a fixed obstacle in the stream. 

Many attempts have been made to solve this equation analytically, but none has 
met with complete success. Most of the progress has been through some linearizing 
process, such as writing g as a power series in M?, the non-linear terms on the right 
hand side of (1.1) then being of higher order in M? than the linear terms on the left. 
Again, a transformation to new independent variables, the magnitude and direction 
of the local fluid velocity, gives a linear equation for gy. In this paper, use is made of 
an expansion of ¢ as a power series in r—! and it is shown that the equation for each of 
the coefficients in the series is linear. Some light is thrown on the form of the expres- 
sion for g which may be expected in the two cases of zero and finite circulation about 
the obstacle. 

2. Case of zero circulation. For the case of flow past an obstacle about which there 
is no circulation we make the assumption that ¢g can be written in the following form: 


= rcos + + r-*f3(0) + ---. (2.1) 


It will be seen later that the term independent of 6 reduces to a constant when the 
circulation is zero, and thus may be ignored. Whether or not this expansion for ¢ is 
valid will be apparent from the subsequent solutions which are obtained for the func- 
tions f,(@). Invalidity of the assumption (2.1) may be made manifest either as a lack 
of convergence of the power series in r~', or as a lack of periodicity in 6 of any of the 
fn(0). In the former case, the series could presumably be made convergent by choosing 
r large enough, so that any conclusions would apply only to regions far from the 
origin. If, as seems probable, there is no specially favoured region, we may regard the 
expansion (2.1) as valid generally provided that the resulting equations for the f,,(8) 


give periodic solutions. 


* Received May 15, 1944. 
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The expression for g? corresponding to (2.1) is 
g 2f, cos — siné) +--- 

+ 2nf, cos — 2ff (2.2) 
where dashes denote differentiation with respect to 8, and, in the general coefficient, 
only terms involving f, have been retained. Substituting both (2.1) and (2.2) in the 
equation (1.1) and equating coefficients of r~"~? we obtain, after some reduction, 

— M? sin? — + 1)M? sin 20 + f,[n? + nM? — (n? + 2n)M? cos? 6] 
M°E,(fn—2, fi, 6), (2.3) 


where E, is an undetermined function. 
The equation (2.3) may be normalized by writing 


h,(8) 
= (1 — M2? sin? (2.4) 


which gives for h,(@) 
(n? — 1)(1 — M?) 
(1 — M? sin? 6)? 


ha’ + at + = M*E,(1 — M? sin? (2.5) 


The homogeneous equation may also be written in the form 


ht + in(1 + uo + >> 2u, cos 200) = 0, (2.6) 
p=1 


where 


2 


up = (n? — 1)(1 — M*)(— 1)?>> (2.7) 


m=p (m — p)\(m + p)! 
The expression for u, is convergent provided M?<1 and }-u, is absolutely conver- 
gent for the same condition. The equation (2.6) will be recognized as being identical 
with Hill's equation which arises in a determination of the motion of the lunar 
perigee.’ As is already known,? Hill’s equation may be solved by the assumption of a 
series expansion of h,; for instance, we may write 


= (2m) !(m + 1) 


hn = (2.8) 
and determine the constants 6,, and wu, which governs the periodicity of the solution. 
However, the method is tedious and we shall find it more expedient to derive a solu- 
tion in powers of M? from the equation (2.3). 
We may also note that if only terms of order M°® and M? are retained in (2.6), 
the homogeneous equation for hk, degenerates to Mathieu's equation, viz., 


hi + h,[n® — M?(n? — 1) cos 26] = 0. (2.9) 


It is known that Mathieu’s equation yields a periodic solution under certain condi- 


1G. W. Hill, Acta. Math., 8, 1-36 (1886). 
2 See for example, Whittaker and Watson, Modern analysis, Camb. Univ. Press, Cambridge, 1940, 


p. 413. 
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tions; a proof that equation (2.9) gives a periodic solution (with consistent approxima- 
tion) is demonstrated in the Appendix. 
3. Periodicity of f,. Consider first the homogeneous equation for f,, viz., 


’(1—M? sin? 6) — (n+1)M? sin 20+, [n?-+-nM?—(n?+2n)M? cos? 6]=0. (3.1) 


In order to show that this equation gives a periodic solution for f,, we write the solu- 
tion in the form 


= + MF + + (3.2) 
Substituting in (3.1) and equating coefficients of powers of M? we obtain a series of 
equations, each of the same type, for F,,0, - . The typical equation is 


nti WF = sin? + + 1) sin 20 — Fajm[m — (n? + 2n) cos?6] (3.3) 


the solution of which is periodic unless a term cos n@ appears on the right-hand side. 
Suppose that F,,,, consists of the sum of a number of cosine terms; if, for instance, 


= C cos (p60 + 3) 
then (3.3) becomes 
+ = + n)(p + + 2) cos [(p + 2)6 + 4] 
+ 4(n? — p*) cos (p0 + 8) +3(n — p)(n — p + 2) cos [(p — 2)0+4]}. (3.4) 
F,,, m41 Will thus be unperiodic only if p=n—2; the cases p=m and p=n-+2 do not 
give finite terms in cos 76 on the right of (3.4). Now the value of Fo is 
A,, cos (nO + €,) (An, €n arbitrary). 


F,,, therefore contains terms in cos (n+2)@ and cos 76. Similarly, F,,2 will contain 
terms in cos (n+4)6, cos (n+2)@ and cos n@. Generally, we see that a term cos (m —2)6 
cannot arise in any of the F,,m. Hence each of the F,m41 will be periodic, and the solu- _ 
tion of (3.1) is periodic. 

Unfortunately, a proof that the particular solution of the equation (2.3) for f, is 
periodic is not as simple, since the function EZ, is a very complicated one. For both 
n=1 and »=2, E,=0 so that the complete solutions of f; and fz are periodic. 

However we can show generally that to the order of M?, at least, f, is periodic. 
To this order the denominator of the right hand side of equation (1.1) may be treated 
as unity. It is also possible to write, in the numerator of this expression, 


fn = Ancos (nO + €,) (3.5) 


where A, and ¢, are constants independent of M*. Considering only the terms of EZ, 
which can lead to a term in cos 6, we find that 


n—2 


cos p(n — p)(n — p — Ifpfr—p-1 — 2 cos (n — fn—p-1, 


n—2 


— 2sin0>> p(n — + 2 sin >> (n — p)(n — p — fn—p-1 


+ unimportant terms. (3.6) 


By the use of (3.5), this becomes 


ot 
4 
— 
ou] 1 
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n—2 
E,= —2 cos >> p(n—p)(n—p—1)A pAn—p-1 Cos 
1 


n—2 


showing that no terms in cos n@ occur in E,. Hence f, is periodic at least as far as the 
order of 

By carrying out in detail the iteration process for a solution of (2.3) in powers 
of M? for fs and f,, and also for f, for the case of a doubly symmetrical body (for which 
fo=fs= +++ =0), the author has been able to show that non-periodicity does not 
arise when terms of the order of M* are considered, nor does it arise when the order 
of M® in f; is considered. In addition, the solution of equation (1.1) appropriate to 
flow past a circular cylinder has been found* to be of the form (2.1), at least as far 
as the order of M®. 

Hence, although not yet proved, it seems probable that the particular solution 
to the equation (2.3) for f, is periodic. If this be the case, then the complete solution 
for f, is periodic and the assumption of a series expansion (2.1) for g may be regarded 
as valid for the case of zero circulation. 

4. Numerical values of the constants in f;. In the case m =1, the normalized equa- 
tion (2.5) becomes 


hi’ +h=0 (4.1) 
which gives 
A cos (@+a 

(4.2) 
1 — M?* sin’ 0 


where A and @ are arbitrary constants, to be determined by the inner boundary layer 

condition. This appears to be the only case for which the complete solution for f, is 
obtainable, and has also been put foreard by Imai.‘ 

In general, A and @ are functions of M* and we may write, for small values of M?, 

A = a) + a,M* + a.M*+-:-, (4.3) 

Qa ao + (4.4) 


Expanding f; in powers of M? we then have 
fi=ao cos (0-+-a) +M?[ao cos sin? 6+-a; cos(8-+a) — aoa sin(@+a») |+ (4.5) 


The value of f; correct to the order M? is known for several different obstacle shapes 
from the work of previous investigators. For instance, Kaplan® has determined the 
first approximation to compressible flow past an elliptic cylinder at an angle of inci- 
dence 8 to the uniform stream. Comparing his result with (4.5) we find, after some 
reduction, 


3 I. Imai, On the flow of a compressible fluid past a circular cylinder, 11, Proc. Phys.-Math. Soc. Japan, 
23, 180 (1941). 

I. Imai, loc. cit. 

°C. Kaplan, On the use of residue theory for treating the subsonic flow of a compressible fluid, N.A.C.A. 
Technical Report 728, 1942. 
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c? sin 28 
a=1- $c? cos 28 + (4.6) (4.7) 
1 
1 — [(— pect — 4 — 16c~*) + (3c? + 20c-*) cos 28 + (qgc* — 3)cos 48] 
0 


= [ae + 14c-? + 32c-*) + (— — 5 — 40c~*)cos 28 


1 4 

+ log = 
a 

+ (3c? + 6c-*)cos 48], (4.8) 


4— 


1 4+? 
a, = sin 28] + + —— (1 — 3c? cos 28) + $(4c~* — 3c?) log . (4.9) 
2a? 4—¢ 
where c is the distance between the centre and focus of the ellipse, and the unit of 
length is the radius of the circle into which the ellipse transforms (equal to the mean 


of the major and minor axes). At zero incidence, these expressions become 


2t 


ao = (4.10) 
4+ c? 
a; = — 1)(3 — + — 1)4 log 
“1 — 2¢ 2t4 
= log (4.11) 
(1 — #)? (1 — #)(1 — 
a =a, = 0, (4.12) 
Ke) 
08 
0.6 
45:4, 
04 7 
02 
- 
0 0.2 0.4 06 0.8 1.0 
THICKNESS RATIO t 
IFLAT PLATE CIRCLE 


Fic. 1. do, a, for ellipses (0° incidence). 


where ¢ is the ratio of thickness to chord of the ellipse. The expressions (4.10) and 
(4.11) are shown graphically in Fig. 1. For ellipses of thickness ratio less than 0-2, 
a, is given accurately by a;=4. 


| 
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Using results obtained by Tomotika and Umemoto,‘ it is also possible to derive 
numerical values for a, a; corresponding to the flow past a symmetrical Joukowski 
aerofoil of arbitrary thickness at zero incidence. After considerable reduction, we 
obtain from a comparison of these authors’ formulae with (4.5) 


a = 3(1 — k)(3 + &), (4.13) 
a, = — — 3k + 2k? + — + + &)® log (1 — 
+ — 12k-5 + 15k-* + — — + 207 — 170k 
— 61k? + 36k* + 10k* — 4k), (4.14) 
a = a, = 0, (4.15) 
where 
1—e 
k= (4. 16) 


and ¢ is the parameter controlling thickness-to-chord ratio in the Joukowski trans- 
formation as shown in the table. The unit of length is again the radius of the circle 


TABLE—Relation between « and ¢. 


e=0 0-03 0-05 0:07 0-10 0-15 0-20 0-30 0-50 
t=0 0-0378 0-0618 0-0849 0-1179 00-1687 0-2150 0-2958 0-4210 1 


into which the aerofoil transforms. The values of ao and a; are also shown in Fig. 2. 
Most practical aerofoil shapes have a thickness ratio less than 0-2 and it is evident 
that ao and a; may be estimated with fair accuracy for these shapes by assuming a 
linear dependence on ¢. Using the slope appropriate to =0, we have 


a= 1.5382, a, = 0.769% = (4.17) 


From the closeness of the values of a» and a; for an ellipse and a Joukowski aero- 
foil of the same thickness ratio, it seems reasonable to suppose that a» and a; do not 
vary greatly with change of shape of the obstacle. These results may then be of use 
in estimating compressible flow past obstacles of different shapes at points far from 
the origin. 

The same procedure is also possible with the higher coefficients in the series for ¢. 
For instance, the equation for f2 is 


fi (1 — sin? 6) — M? sin 20 + f2(4 + 2M? — 8M? cos? 6) =0 (4.18) 
of which the solution, in powers of M*, as far as M’, is 
fa = bo cos (20 + Bo) + M*[b: cos (20 + B:) — $bo cos (40 + Bo)] +--+, (4.19) 


where Do, 5;, 80, 6: are arbitrary and independent of M*. These constants could also be 
determined from a comparison with known first-order solutions of the compressible 


6S. Tomotika and H. Umemoto, On the subsonic flow of a compressible fluid past a symmetrical 
Joukowski aerofoil, Tokyo Imp. Univ. Aero. Res. Inst., Report 205, 1941. 

The expression for the velocity potential for compressible flow given by Tomotika and Umemoto is 
slightly in error. 

In the original symbols the coefficient of cos (25+) in the expression for g (Eq. 122) is given as 
(h®/k)(2 + d» + + 23 — 8d2 log Ax) but should read as (h#/k)(2 + + — 8A2 log As). 
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flow past particular obstacles. We may note that for bodies symmetrical about axes 
parallel and perpendicular to the uniform stream, 


bo = = --- = (4.20) 
For these doubly-symmetrical bodies, in fact, 
(4.21) 
1.0 
0.8 


"4 


0.2 
02 O04 O68 O08 10 
THICKNESS RATIO t 
FLAT PLATE CIRCLE 


Fic. 2. ao, a; for Joukowski aerofoils (0° incidence). 


5. Case of finite circulation. In the case of flow past a body about which a finite 
circulation exists, the assumption of a power series expansion for g becomes 


= + fo(6) + + +--+, (5.1) 


where fo is not periodic, but fi, f2, - - - must be so. Once again we can take the cri- 
terion of validity of this assumption to be the periodicity of the resulting solutions 
for fi, fe, - - +. Substituting (5.1) in the original equation (1.1) and equating coeffi- 
cients of r—*~?, the equation for f, is obtained as 
fi (1 — M? sin? 0) — f,) (n + 1)M? sin 20 + f,[n? + nM? — (n? + 2n)M? cos? 6] 

(fr—1) Sn—2; fo 6), (5.2) 
where 


E, = E, — 2(n + 1) cos 6fof.—1 — 2 sin 0fo f.-1 — 2 sin Ofofa—s 


p=0 


— (y — — + 
+ — fa-p-2] (n> 1). (5.3) 


‘ 
| 
‘ 
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As before, the appearance of a term cos ”@ in E,! at any stage of the solution in powers 
of M? will indicate that f, is non-periodic and hence that the assumption (5.1) is in- 
valid. We shall in fact show that in this case, non-periodic solutions do arise. However, 
a more appropriate form for the function ¢g presents itself, and this is shown to lead 
to a solution which is self-consistent. 

Confining attention to the order M?, we may write 


fn = Fao + M*F,,1. (5.4) 
F,,,0, of course, is the incompressible solution 
= An cos (nO + €n). (5.5) 


Substituting (5.4) and (5.5) in (5.2) and equating coefficients of powers of M? we ob- 
tain the equation for F,,, as 


Fy + = Agn(n + 1) cos + + + An—a{2n(m — sin (nO + 
+ 2(n — sin [(m — 2)0 + en-1] 
+ 2(n — 1)fé" sin @ sin [(m — 1)6 + 


n—2 
— — p — 2) sin [(m — 2p — 2)0 + €n—p-2 — } 
p=0 
+ contribution from E,. (5.6) 


It will be shown later that f¢’ is of order M* and may thus be neglected here. As shown 
previously E, does not give rise to a term cos 6 or sin v6 in this equation. Hence there 
is one term on the right of (5.6) which gives a non-periodic particular solution for Fp.,1; 
this is 
2n(n — sin (nO + €n_1) (5.7) 


and the power series expansion (5.1) for g is not valid. 
A form for ¢ which will avoid an inconsistent solution is suggested by the relation 


2 0G 
[log r-G(r, 6)] = — + log r-V°G, (S.8) 
r or 


where G is an arbitrary function of r and 0, and V? is the Laplacian operator. If we 
choose G to be a solution of Laplace’s equation, in particular, 


G = B,:r™ cos + (5.9) 
where B, and e€, are arbitrary constants, then 
log r) = — cos + ). (5.10) 
If now, in place of (5.1) we write 
fot + log r-M?B, cos (nb + (5.11) 
n=1 


the equation (5.6) for F,,; is unaffected except for the addition to the right-hand side 


of a term 
2nB, cos (nO + ). 
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The particular solution for F,,, will thus be periodic if we choose 
(n> 1). (5.12) 
The case of »=1 warrants special attention since we cannot write Fy_1,0 


=A,-1 cos [(n—1)0+¢,-1] when n=1. When 1, the “dangerous” term on the right 
of (5.6) is — 2f/? cos 6 so that we must choose B,=fj?, e{ =0. Hence, 


>> r-"B, cos (nO + €,) = fo > 


ami r =" 
CoS [A, cos (nO + €n) 
= + fo {00s 0 | 
1 700 
TA, cos * 
+ sin (5.13) 
Or 
But 
r cos 0 + + >. cos (nO + €n) = go, (5.14) 
n=1 


where ¢ is the velocity potential for incompressible flow and 27k (=2zf¢, neglecting 
terms in M/”) is the corresponding circulation. Thus 
ag ag 
>> r-*B, cos (nb + = ko (cos + sin 6 = (5.15) 
1 700 or oy 
where y is the space ordinate at right angles to the direction of the free stream. 
We have then, that when a finite circulation exists, the expression for ¢, correct 
to the order M? at least, is ; 


g=rcos0+ fat log (5.16) 
n=1 y 
and this expression is confirmed by the known results for a circular cylinder.? We may 
also note, as a matter of convenience, that the expression for f, in (5.16) may be ob- 
tained by formally ignoring the log term when substituting in (1.1), and rejecting any 
non-periodic solutions. 
6. Solutions for f) and f;. Putting »=0 in (5.2) and noting that Ey =0, we have 


M? sin 26 : 
1 — M?sin?0°°’ 


0 


integration of which gives 
K 


j= 6.2 

fo 1 — M? sin? 0 
where « is a constant which may depend on M?. This equation shows that the circula- 
tion about the body in compressible flow is 


f fo (1 M?)}!2 (6 ) 


7S. Tomotika and H. Umemoto, loc. cit., Appendix. 
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Substituting for « from (6.3) in (6.2) gives f¢ in terms of the circulation, and an 
identical equation for the transverse velocity at large r has been found by Glauert*® 
using “small perturbation” theory. 

The contribution to the velocity potential is found by integration of (6.2) to be 


fo = — M*)—/? arctan [(1 — M%)"/? tan 6] 


or, in terms of the circulation C, 
fo= arctan [(1 — M?)"/? tan 6]. (6.4) 


Just as the magnitude of the circulation about a body in incompressible flow is arbi- 
trary, so in compressible flow, the relation between circulation and Mach number is 
arbitrary. In certain cases an independent criterion for this relation is available, e.g., 
for sharp-tailed bodies the circulation should be such as to prohibit an infinite velocity 
at the tail-point for all Mach numbers. For blunt bodies, a convenient criterion is 
that the circulation shall be the same at all Mach numbers. In this case, we may treat 
C as a constant and the expansion of fo in powers of M? is 


fo= [@ — 4M? sin 20 + } sin 26 + sin 40) + --- ]. (6.5) 
T 


(Actually this expression for fo has not been proved valid for powers of M? beyond 

the first. However it seems probable that equation (6.1) for fo will hold whatever the 

form of the remainder of the expression for g. In the case of a circular cylinder with 

circulation, it is not hard to show that (6.5) is correct to the order M‘* at least.) 
Putting m =1 in equation (5.2) ,we have 


i’(1 — M? sin? 6) — 2f{ M? sin 26 + fi(1 + M* — 3M? cos*?@)=0, (6.6) 
the only finite terms in Ey being terms which would give rise to non-periodic particu- 
lar solutions. As shown in section 2, the solution is 

A cos (6 + a) 
1 — M’ sin’? a 
where we are now justified in using this expression to the order M? only; i.e., 


fi = a cos (0 + ao) 
+ M?[ao cos (@ + a) sin? @ + a; cos (0 + ao) — aoa; sin (6 + ap) ], (6.8) 


the constants having the same meaning as in section 2 but different values owing to 
the finite circulation. 

As in the previous work, it is possible to derive the values of ao, a1, ao and a, by 
comparison of (6.8) with the known results for certain obstacles. Results are available® 
for the flow of a compressible fluid about a circular cylinder of unit radius, with con- 
stant circulation C. Comparison with (6.8) shows that 


a= 1, a= §+ (6.9) 


* H. Glauert, The effect of compressibility on the lift on an aerofoil, Proc. Roy. Soc., A 118, 113 (1928). 
* S. Tomotika and H. Umemoto, loc. cit., Appendix. 
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so that the value of f; is not greatly altered by the presence of a small circulation. 

We may also note that for a body, symmetrical about axes parallel and perpendic- 
ular to the uniform stream 

fi (6) fi(x = 6), 

giving a =a,=0 and a corresponding result holds for fs, fs, - - - . In the case of even 
values of m, the angle constants in f, are equal to t. 

7. Appendix; Proof that the equation y’+y[n?— M*(n?—1) cos 20]=0 gives a 
periodic solution. The general form of Mathieu's equation!® may be taken as 


y” + y(a + 169 cos 26) = 0 (7.1) 


the solution of which is periodic only if a and g are suitably related. The required 
relation may be written 
a=n?+ ag + (7.2) 


since a must reduce to the square of an integer when g vanishes. Now, for the particu- 
lar case of Mathieu’s equation with which we are concerned, 


a= n’, 16g = — M*(n? — 1), (7.3) 


and (7.2) will be satisfied to the order of M? if a;=0, or, when n=1, by any value 
of a;. Thus we need to show that a,=0 for the general Mathieu equation. 
This may be done by assuming a solution of (7.1) of the form 


y = cos nO + ga;(0) + +---, (7.4) 


where a;(@), a2(@), - - - are periodic functions of 0, independent of g. (The proof for 
solutions which reduce to sin 8 when g=0 is identical.) Substituting (7.2) and (7.4) 
in (7.1) and equating coefficients of g we obtain 


al + na; + cos n6(a, + 16 cos 26) = 0, 
i.e., 
ay + n’ay = 8 cos (n + 2)0 + a, cos nb + 8 cos (n — 2)0 
giving a:=0, (#1) since a is periodic. For the particular case m=1 we have instead 
a,= —8, but in this case we no longer require a; to be zero. 


Further coefficients a2, a3, - - - and the functions a;, a2, - - + may be determined 
by equating higher powers of g; the general Mathieu function may thus be constructed 


if desired. 


10 Whittaker and Watson, loc. cit., p. 409. 
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THE RESONATOR ACTION THEOREM* 


BY 


W. R. MACLEAN 
Polytechnic Institute of Brooklyn 


1. The problem. In the kinetic theory of gases there are a series of theorems relat- 
ing to the “adiabatic invariants” of mechanical systems.! One such theorem, which 
is a special case? of a more general theorem of Boltzmann,’ is contained in the state- 
ment: 

In a frictionless machine having a periodic motion, the kinetic action, 1.e. the product 
of mean kinetic energy and period, is invariant against an adiabatic deformation. 

The word “adiabatic,” relating to a machine, seems a strange one to use, but the 
meaning here bears an analogy to thermodynamics. If a machine moving in periodic 
motion has additional coordinates which are normally held fixed, but can, at will, be 
varied by outside influences, a variation of these will, in general, change the period. 
An example would be a violin string which is tightened while vibrating. Adiabatic 
means that the deformation is carried out with perfect uniformity and infinite slow- 
ness so that the energy so introduced has only d-c components and hence none near 
the frequency of oscillation. 

If the machine is linear, that is, can be represented by linear differential equations 
(which will have constant coefficients when the machine is not being deformed), the 
theorem can be simplified to: In a linear frictionless machine having a periodic motion, 
the action, i.e. the product of total energy and period, is invariant against an adiabatic 
deformation. The simplification follows from the fact that in this type of system the 
mean potential and kinetic energies are equal. 

The proof of this theorem, as carried out in Jeans and elsewhere, although quite 
general for a machine is still limited to a finite number of degrees of freedom. For the 
case of an electromagnetic resonator (where an infinity of degrees of freedom is in- 
volved) the proof for a cubical box has been indicated by Ehrenfest‘ in a paper on 
quantum problems. The theorem is, however, usually accepted as true in general,® 
but as far as the author is aware, a proof for a resonator of any shape and any dis- 
tribution of dielectric constant as is needed for application to radio problems is still 
lacking. 

It is the purpose of this paper to carry through such a proof for the specific case 
of a general lossless electro-magnetic resonator. Such a proof is considered worth while 
for engineers, since the theorem has certain uses in the design of equipment. For in- 
stance, if one is interested in the frequency change caused by a slight deformation 


* Received June 12, 1944. 

1 For quick reference see for instance J. H. Jeans, The dynamical theory of gases, 4th ed. Cambridge 
University Press, 1925, p. 410 ff. 

2 P. Ehrenfest, Proceedings of the Amsterdam Academy, 16, 591 (1914). 

L. Bolzmann, Verlesungen aber Mechanik, II. J. A. Barth, Leipzig, 1904, §48. 

4 P. Ehrenfest Annalen der Physik, 36, 91 (1911). 

5 For instance, in L. Brillouin, La Statistique Quantique, Les Presses Universitaires de France, Paris, 


1930, p. 39. 
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of an electro-magnetic structure and if this proves difficult to evaluate directly, one 
can instead try to determine the work done against the electro-magnetic forces, which 
by this theorem will give the desired result. 

2. Preliminary observations. The basis of the proof is formed, of course, by 
Maxwell’s equations. For an uncharged region without permanent magnetism, con- 
ductivity, or driving forces, they become.” 

Primary Equation: curl H = D’+G. But here, G = 0. 
Secondary Equation: curl E = — B’. 

Magnetic Equation: div B= 0. 

Electric Equation: div D =0. 

Dielectric Equation: D = «€E. 

Core Equation: B = wH. 

For later purposes, it is convenient to introduce a vector potential in a manner 
somewhat different from the usual one. We suppose that there exists a field satisfying 
the above equations, and then by virtue of the electric equation, find a U such that: 


D = curl U. (1) 


This designation still leaves U in large part arbitrary. From the Primary Equation, 


one can write: 
H — U’ = — grad y, 


where y is a suitable scalar. One now chooses 
div U’ = div H. 

U is still arbitrary to the extent of the gradient of a harmonic (in the sense of poten- 
tial theory) function. But from the last two: 

Ay = 0; 
hence y is a harmonic function. So we can still adjust U so that: 

H=U"’. (2) 

The Secondary Equation finally subjects U to a differential equation: 
curl [(1/e) curl U] = — up”. (3) 


Since we are not restricting ourselves to homogeneous media, this cannot be reduced 
to the usual Wave Equation. However, a step in the proof can be taken by the usual 
device of setting 


U(x, y, 2, 4) = S(x, y, 2) f(d) 
which leads to 
f(t) curl [(1/e) curl S(x, y, z)] = — uS(x, y, 


Since this is identically true in x, y, 2, t, it can only hold if f and f’’ are proportional. 
In particular 
6 The author feels that Maxwell’s equations should have names if for no other reason than to facilitate 


reference to them. “Primary” and “secondary” remind one of the phenomena in an unloaded transformer. 
7 Primes indicate derivatives with respect to time. 
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wf, 


the solution of which is: 

f = a cos wt + sin ot. 
Hence, individual solutions of (3) are sinusoidal in time. Since (3) is linear, solutions 
are made up by a superposition of such particular solutions. Substituting the last into 
(3) leads to: 
curl [(1/e) curl S] — wuS = 0. 

This will have solutions satisfying the boundary conditions (E normal and B tan- 
gent to the metal surfaces) only for certain values of w, the so-called eigenfrequencies 
of the resonator. In general, these eigenfrequencies will not be harmonically related. 
Since the theorem requires that the resonator have a period both before and after 
the deformation it is necessary to restrict our considerations to the case where only 
one mode is excited at a time. We must have only one mode, i.e., one frequency also 
after the deformation. This precaution is made necessary by the existence of the spe- 
cial case of degenerate modes. 

If a given resonator had a set of eigenfrequencies, one would expect after a very 
small deformation that it would then have a different set of eigenfrequencies which 
could be put in one to one correspondence with the old ones, and wherein correspond- 
ing ones would differ by only a very small amount. In particular, if the deformation 
were made small enough, these differences could be made small compared with the 
spacing of the eigenfrequencies of the original or final resonator. 

This expectation is correct with one exception. Suppose that a particular mode Ma, 
having radian frequency w, were excited. Some field quantity, say the potential U 
would be a particular function of position :* 


U = U.(x, 2). 


In general, this is the only function which can satisfy the differential equation and 
the boundary conditions. As such it is called an eigenfunction of the problem belong- 
ing to the eigenfrequency, w,. As an exceptional case, however, it can occur that more 
than one eigenfunction, i.e., more than one mode exists for the same frequency. Then 
this particular frequency is called degenerate. The number of linearly independent 
modes corresponding to one frequency is called the degree of degeneracy. This phe- 
nomenon is mathematically similar to the Stark and Zeeman effects. 

A degeneracy depends on the existence of some type of symmetry in the resonator. 
For instance, in a cylindrical resonator, the TE: mode® is degenerate. 

A single TE; has a transverse axis of orientation, say the diameter along which E 
is radial. One can, however, superimpose thereon another TE, orientated at right 
angles to the first and having perhaps a time phase angle against the first. They are 
obviously linearly independent and have the same frequency. If now such a cylindri- 
cal resonator were deformed in such a manner as to destroy the symmetry, the de- 
generacy would in general be broken and the frequencies split. As a result, the condi- 
tion that the resonator have a period after the deformation would not be met. 


8 The bar over a letter indicates complex vector phasor. This is the usual notation for simple harmonic 


fields. 
® This designation of modes is now common. See, for example, Sarbacher and Edson, Hyper and ultra- 


high frequency engineering, J. Wiley and Sons, 1943, p. 387. 
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In this special case, the theorem applies to each component separately. The proper 
component in the degenerate state can always be chosen as follows: one deforms the 
resonator splitting the frequencies, then quenches all but one frequency, then relaxes 
the deformation. Repetition of this deformation will then yield on!y one frequency. 

For purposes of the theorem, it is necessary to demonstrate that the mean electric 
and magnetic energies of a resonator oscillating in one mode are equal. In complex 
amplitude notation one has for the magnetic energy:!° 


2(Wau) = f dv H*-B. 
V 


Applying the Secondary Equation and the vector identity for the divergence of a cross 
product: 
div (A X B) = B-curlA — A-curlB (4) 


one has 
2(Wmu) = f dv H*-j(curl E)/w = div (E X H*) + f (jdv/w)E - curl H*. 
Vv V 


The first integral on the right transforms to the surface by Gauss’ theorem and 
vanishes by the rule of exchanging dot and cross in the triple product since E is nor- 
mal to the surface. In the second integrand we use the primary equation and have 


2(Wau) = = 2(Wz). 


If one solves for the field, say for U, in a resonator oscillating at one frequency, 
it is conceivably possible for U to be a complete complex vector, i.e., have non-colinear 
real and imaginary parts. If this is so, however, then the mode is certainly degenerate, 
since the real and imaginary parts of U each satisfy the differential equations. The 
fields derived from each part of U individually satisfy Maxwell’s equations and the 
boundary equations. So as long as we restrict ourselves to non-degenerate modes it 
suffices to consider only field vectors which swing on a line. Moreover, since only pure 
real or imaginary U is needed, U is everywhere in phase, as are E and H. The latter 
two are, however, in exact time quadrature with each other as is seen from Maxwell’s 
Equations. 

To compute the work done in deforming the walls of the resonator or moving the 
pieces of material within it, one needs the expression for the density of force exerted. 
For the material, this is found by taking the divergence of the Maxwell tension 


tensor," 7: 
Ti; = ED; — + HB; — 35;H-B. (S) 


The density of force, k, is then found from 


3 
ky = AT ;;/dx; 


j=1 


10 The marks < > designate time average over a cycle. Asterisk means complex conjugate. The com- 


plex amplitudes are rms. 
1 §;; is the Kronecker delta: equals unity if the indices are alike, otherwise zero. Note we are now 


working with instantaneous values. 
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which can be evaluated to give 


d 
k = — grad — grad + XB) 


for an uncharged, passive, and lossless region. 

Since the deformations to be made are to be adiabatic, i.e., carried out uniformly 
with extreme slowness, the work done can be obtained from the mean value of the 
force. Averaging the last over a cycle, one has: 


(k) = — 3(B*) grad « — }(H?) grad u. (6) 


This is the density of force exerted by the field on the material. 

To find the force exerted by the field on the walls of the resonator, pick a certain 
point on the walls and there make a coordinate transformation bringing the x axis 
coincident with the inward normal. E has then only an x component. At the same time 
make the y axis coincide with the direction of the magnetic field which lies in the 
boundary. The tension then takes on the form 


WE 0 0 —wu 0 0 
Ti; =| 0 -—-we O J+ 0 we oO |, 
0 0 — Wr 0 0 — WM 


where the w’s are the electric and magnetic energy densities. Since the vector element 
of surface has only an (positive) x component (for inward normal) the surface density 
of force, f, becomes 
f = (wz — wy), 0, 0. 
Since the x component is directed inwardly, we can say that the field exerts either 
pure pressure or pure suction on the walls. In particular, if we designate as the density 


of the Lagrangian function: 
l= wy — WE; (7) 


the field exerts an (algebraic) pressure, p, on the walls, given by: 
p=l. (8) 


For slow uniform (adiabatic) movements of the walls, the work done can be 
found from the mean pressure 


1 
= — f lat = (9) 


where Tf is the period of oscillation. 

3. The resonator theorem. Jn a lossless electromagnetic resonator, the action of each 
mode, 1.e., the product of total energy and period, is invariant against an adiabatic de- 
formation. In case a symmetry creating a degenerate frequency is involved, the separa- 
tion of modes is to be carried out as described above. 

Since the mean electric and magnetic total energies are equal, the integral of / over 
the volume V and the period 7 vanishes: 


h [a fii = 0. (10) 
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Since this is true for any resonator, it is also true after we have deformed our resonator 
slightly. Hence, the variation of h, i.e., the difference between a new and an old hk, 
is zero. 

To see most clearly what is meant, it is best to imagine two editions of the original 
resonator at hand, one of which is varied. After a slight variation, they beat ever so 
slowly with each other. Since throughout each, U is in phase, one waits for a cycle 
when the U in both resonators pass through zero together. This instant we take as 
t=0. Since they are almost alike, the two U’s will differ by only a small amount over 
the following cycle. 

The variation in h, 6h, which is zero, is brought about by four effects: (a) the 
change in period, (b) the change in volume, (c), the change in the field, (d) the change 
in dielectric constant and permeability due to movement of the materials.” These 


four increments together lead to 


5h = br[L], + rf f at = 0, 
bv 0 


where L is the volume integral of /, i.e. Z is the Lagrangian function. Since L is peri- 
odic, the initial rather than final value may be used in the first term. The second term 
is integrated over 5V which consists of a thin shell adjoining the envelope, hence the 
value of </> to be used is that obtaining on the closest surface point. One sees, there- 
fore, that the second term is proportional to the work done by the field on the walls. 
Hence, in terms of the field energy, W, it is 7 times the negative increment of W due 
to the movement of conductors. So the last becomes: 


r 1 1 
0 = — 7(6W), +f at f dos (— ae ~ =p’). 
0 2 € 


In the expression for /, the magnetic and electric energy densities are written in 
inverse form, so to speak, for reasons of future convenience. Carrying out the indi- 


cated variation in the last term above one gets 


1 
0 = — +f at — E-5D)+ —f at f + E%ée). 
0 V 0 V 


Equation (6) gives an expression for the mean value of the force density on the 
material. If one multiplies (6) by 6x, i.e. by the displacement of the material, and in- 
tegrates over the volume, one has the work done by the field on the material. This is 
proportional to the last term above, the latter being just 7 times the negative incre- 
ment of field energy due to movement of the material. Also by (2) the quantity H 
in the first integral can be expressed in terms of the potential, U, leading to: 


0 = — 7(5W). +f at -f at — 
0 V 0 V 


The variations of W due to conductor and material movement combine to form 


22 The question of surfaces of discontinuity on material boundaries is ignored. This is permissible, 
since they may be replaced by large though finite values of the gradient of the material constants. It is 
merely a matter of mathematical convenience which method is adopted. 
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the total variation of field energy. In the first integral we replace the integrand by 
the time derivative of both factors less the correction term and integrate the former 
with respect to time. We also express D in terms of U by (1), giving 


0 = Loir — ow + f f at + B-cur (11) 
Vv 0 


Since B has period 1, it can be taken outside the time limits in the first integrand, 
which leaves 6U between those limits. However, 5U=0 at t=0 by our choice of time 
origin. At t=7, the old U is also zero, hence 6U at that time is simply the new U, 
say U, at =r. The x component of the latter would be given, for example, by 


U.(t) = A sin (w + dw)t. 
The amplitude can be expressed in terms of the initial derivative, and if we also 
substitute we have: 
1 
6U,|, = U(r) = ——— (0) sin dw r. 
[sU.] (r) (0) T 


For small variations, the sin may be dropped, and if 5w is expressed in terms of 6r, 
one has to first order: 


[65U.], = — UL = — [UZ 


In the last step it is permissible to equate the old to the new since we have already 


obtained the variation. 
Substituting this result into (11) one has: 


0 = Lor — — ir av[B-U" -f at + B-curl 5U). 
Vv 0 


Since U’=H, the first integral is the negative of twice the magnetic energy at 
t=0 times the variation in period. This combines with the first term to give — Wodr. 
But since the total energy does not vary in time, one has by changing sign in the last: 


0 = Wér + r6W +I, (12) 


where J is the last integral above which we will now show to be zero. By the use of the 
identity (4) the last term can be integrated by parts which gives 


I = f fa f dw + 
0 8 0 V 


The second integrand is zero by the Secondary Equation, and the first is seen to 
be zero by interchanging dot and cross since E is normal to the surface. 
Using this in (12) we have simply: 


5(7W) = 0 


which is the theorem. 
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II. Analysis of a Given System with the Help of the 
Characteristic Function, Using the 
Direct Method of Analysis* 


BY 
M. HERZBERGER 
Communication No. 961 from the Kodak Research Laboratories 


In previous papers':** the author has developed a direct method of geometri- 
cal optics better adapted to the actual procedures of optical calculation than pre- 
vious methods. In this paper we shall introduce in this direct theory Hamilton's 
characteristic function, and with it we shall find the image of an arbitrary surface in 
a given optical system with rotational symmetry. The author apologizes for changing 
his notations again. In the previous papers a ray was specified by the coordinates 
(x, y) of its intersection point with the object plane and by the optical direction cosines 
£, n, $=Vn2?—#—7?, i.e., the direction cosines multiplied by the respective re- 
fractive indices. In this paper, x, y represent the optical lengths (i.e., lengths multi- 
plied by the refractive index) of the coordinates of the intersection point with the 
coordinate plane, while £, 7, ¢=\/1—#—7? are the direction cosines. The reader 
may prove for himself that the fundamental formulas given in the former papers re- 
main unaffected by this change in coordinates. 

1. Hamilton’s characteristic function and the direct method. In the previous pa- 
pers we introduced the vector b with coordinates x, y, the vector t with coordinates 
€, n, and analogously the vectors b’(x’, y’) and t’(é’, n’) for the image ray, x’, y’ being 
the optical lengths of the coordinates of the intersection point of the image ray with 
the coordinate plane z’=0, and &’, n’ the direction cosines of this ray. 

We found the equation 


b’ = ab + ft, t’ = yb + bt, (1) 
where a, 8, y and 6 depend only on the symmetric functions 
b? = u, bt = 2, t? = w. (2) 


They are connected by one finite equation, 
ad — By _ 1, (3) 


and three differential equations. 
We can find a first integral of these three differential equations by using Hamil- 
ton’s characteristic function. Hamilton has shown that there exists a function V which 


* Received Feb. 5, 1944. Studies in optics: I. General coordinates for optical systems with central or 
axial symmetry appeared in this Quarterly, 2, 196-204 (1944). 

1M. Herzberger, Direct methods in geometrical optics, Trans. Am. Math. Soc. 53, 218-229 (1943). 

2M. Herzberger, A direct image error theory, Quarterly of Applied Mathematics, 1, 69-77 (1943). 

3M. Herzberger, Theory of transversal curves and the connections between the calculus of variations 
and the theory of partial differential equations, Proc. Nat. Acad. Sci. U.S.A., 24, 466-473 (1938). 
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is a characteristic of the optical path (sum of optical distances) between the object 

and image planes. If we know this function V as function of x, y, x’, y’, we find that 
OV OV 

(4) 


Ox Ox’ dy dy 


/ 


These equations can be combined into the vector equation 
dV = t’db’ — tdb, (S) 


where the left side is a total differential. 
Since (5) does not depend upon the choice of variables, we can therefore assume 
x, y, £, » to be independent variables, and replace b and t by their expressions in (1). 


We find then 
dV = (yb + St)(adb + Bdt + bda + tdB) — tdb 
= }[aydu + 2Bydv + Bidw] + y(uda + vdB) + 5(vda + wd). (6) 
The reader may convince himself that the differential equations in the former paper 
are simply the three integrability conditions of (6), and that we obtain them by cal- 
culating Eu», Euw, Eow from (6). Equation (6) is equivalent to 
Vu = fay + y(wau + Bu) + 5(vax + 
Ve = By + + Br) + + wr), (7) 
Vw 386 + + vB) + 5(vaw + wBw). 
We can eliminate one of the four coefficients, a, 8, y, or 5, using (3). Eliminating 7, 
we are led to 
BV, = 5[4a? + (au + + (av + Bw)B.] — (fa + + 8.2), 
BV, = 5[ap (au (av + | (8 + anu + B.2), (8) 
BV» = 5[38? + (au + Bo)aw + (av + Bw)Bw] — (awa + Bw). 
We can eliminate 6 from (8), and thus get two differential equations connecting @ 
and B. If V isa given function of u, v, w, we can thus (theoretically) determine @ and 8, 


and therefore y and 4, by means of this single function V. The differential equations 
in question are 


BV. + fa + + BV, + B+ au + 
ha? + (au + Bo)ay +(av+ Bw)B, a8 + (au + Bo)a, + (av + Bw)B, 
+ + Bud 
1B? + (au + Bo)aw +(av + 


(9) 


or 
VulaB + (au + Bra, + (av + Bw)B.] — + (au + + (av + Bw) Bu] 
— — + + + — aBu)(uw — = 0, (10) 
V.[38? + (au + Bo)aw + (av + Bw)Bw] — VulaB + (au + + (av + 
+ + + Bra + $808, + BwBw + — ewhr)(uw — v*) = 0. 
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Let us now assume that for a given optical system we have calculated V, and there- 
fore a, 8, y, 5; and let us proceed to investigate the image of an arbitrary surface. 

2. Point and diapoint, the diapoint function. In a former paper‘ the author intro- 
duced the definition of a diapoint. Let us consider an object point P and a ray 
originating from it. The plane through the object point and the symmetry axis of 
the system is called the meridian plane. We define the diapoint as the point where 
the image ray intersects the meridian plane, and the diamagnification v as the ratio 
of the distances of a point and diapoint from the axis. The function F describing the 
optical path from the object point to its diapoint shall be called the diacharacteristic. 
To each point of the object space belongs a diacharacteristic, which we shall compute 
later. 

To a given arbitrary point belongs a two-dimensional manifold of diapoints, one 
for each ray, and all these diapoints lie in the meridian plane. The diapoint charac- 
teristic for a given object point thus assumes a two-dimensional manifold of values. 
These may, however, be singular points; for instance, points having a sharp image. 
In this case, all the diapoints fall together, and, according to Fermat’s law, the optical 
path from point to diapoint is constant along all rays. The diapoint characteristic is 
a constant. 

An intermediate case is the one in which all rays intersect a curve in the image 
space. It can be proved that, in this case, the image rays can be split up into a one- 
dimensional manifold of pencils, each pencil converging to a point of this curve in 
such a manner that they form a series of circular cones with the curve tangent to the 
axis of each cone at its apex.’ For each ray in such a cone, the diacharacteristic 
must again have a constant value, and thus the diacharacteristic can assume only a 
one-dimensional manifold of values. Such an image has been called half-symmetric 
by M. Boegehold® and the author. 

Thus, the characteristic describes the image of any object point. Let us investigate 
an object point with the coordinates Xo, yo, 29 (optical lengths) and an arbitrary ray 
through it with the direction cosines £, n, f=V1—(#+n’) =V1—w. Let \ be the 
optical distance along the ray from the coordinate plane to the object point. We find 
that 


A= 2 V1 —~w. (11) 
Let bo be the vector (xo, yo). We then have 
bo + At = b. (12) 
Inserting this into (1), we find that 
b’ = ab+ (B+ad)t, = (6+ yt. (13) 


An arbitrary point on the image ray is given by the vector by with 
bo = b’ + = (a + + (B + ad + + (14) 
We find the diapoint by choosing \’ such that bg and bp are parallel, or that 


4M. Herzberger, A new theory of optical image formation, J. Opt. Soc. Amer. 26, 197-204 (1936). 
5M. Herzberger, Die Gesetze erster Ordnung in optischen Systemen, Zeit. fiir Physik, 45, 86-96 


(1927). 
6H. Boegehold, Raumsymmetrische Abbildung, Zeit. fiir Inst. 56, 98-109 (1938). 
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d’ = (ad — B)/(5 — yd). (15) 
This gives for the diamagnification 
v=atdy = 1/(6 yh). (16) 
The diacharacteristic is now given by 
D=V—-A4+N =V—A+ (adr — B)/(6 — yA), (17) 


where V(u, v, w) is the optical path between the object and. the image coordinate 
planes. D is given in (17) as a function of u, v, w, and z, but for a given object point, 
these are not independent variables. Equation (12) shows that 


Xo — AE = 2, Yo— An= y. (18) 
If we set 
Xo + Yo = Uo, (19) 
then wp is constant for the given object point, and we find that u, v, w are connected 
by the equation 
uo = u + 2d + A*w. (20) 
Equation (20) allows us to eliminate u. Inserting this into (17), we obtain D as a func- 
tion of v and w alone, with uw» and 2 as parameters. This function D(uo, Zo, v, w) is 
the diacharacteristic. We obtain the 
points that are imaged sharply if we 
find the value of uw» and z for which 
the following identities are true: 

D, = 0, Dy» = 0, (21) 
for any value of v and w. If a point 
has a half-symmetric image, we must 
have 


D. Deo @ 0. (22) 


Points uo, 20 for which (22) is fulfilled 
identically have half-symmetric im- 
ages. 

3. Refraction at a sphere. As an 
example let us investigate the image 
formation of a sphere of unit radius. 
The object and image planes will be 
at the center. If s, s’ are unit vectors 
along the ray (Fig. 1), we then have 


Fic. 1. 


a’=aa, s'=ya+is, (23) 
ab = 1, (24) 


We denote the optical path from the object point to the refracting sphere (vec- 
tor a,) by A, and from the refracting sphere to the image reference plane by \’. We 
then have 


with 


a+ As=na,, n'a,+s' =a’. (25) 
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Inserting (25) into (23), we find that 


= — + nn’'yn), a=n'/n+yN, 5 = n/n’ + yh, (26) 
whence 
n? — n'? + nn'yr 
V=A =A+N. (27) 
n? + nn'yr 
Finally, we find from (25) that 
a Xs = n(a, Xs) = n'(a, Xs’) = (a’ X8’). (28) 


Thus, a Xs is an invariant vector. Since a-s =va?=u (cf. (2)), we then have 


na, n'a, = V/n?—u+ (29) 


Equations (25) and (29) now give 


nn 
a= ’ = 0, 
n* + (0 — — Q) (30) 
y= — Q)/mr’, 5 = n/n’ + (Q — — Q)/nn'; 
and for the characteristic function we have 
n* + (0 — — Q) 
Since 
n’? — n? = Q’? — Q?, (31a) 
Eq. (31) can be written in the form 


u+ (Q — + ») 


This is the characteristic function of the sphere. We now form the diapoint charac- 
teristic. Let u be the distance of an object point xo, yo, 20 from the reference plane 
measured along the object ray. Equation (17) gives 


| (33) 
5 — yu (6 — yu)(6 — yd) 
In these equations we substitute for u from the relation 
Uo + po t+ pw. (34) 
If we write 
uo + 2? = p?, (35) 
where p is the distance of the object point from the center, and 
— (36) 
we then find that all our functions depend only on p? and r. We have in particular 
QO? = n? — p? + 7’, = — p? + 7, (37) 


and finally 
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+ (Q — 7)(Q’ + 7) 
F as function of v and w assumes only a one-dimensional manifold of values. Every 


point has a symmetrical image. 
Equation (38) becomes independent of 7 only if 


(38) 


Q=r or Q=-t. (39) 

The first case leads to p=, i.e., the refracting surface itself has a sharp image. The 
second case leads to 

p=—n, (40) 


i.e., to the so-called aplanatic surfaces. In both cases the optical path between the 
object and the image points vanishes. 

A simpler method leads to the diapoint characteristic of a plane refracting surface. 
If z is the distance of the object point from the plane, we find that 


F = 2(1 — w)~/?(n'2/n? — 1). (41) 


For a given object point, F is a function of w alone. Every object point is therefore 
half-symmetric: Only for z=0 can F become independent of w. The plane therefore 
gives a sharp image only of itself (and of the infinite plane, whose points are given 
by 2/f =2/\/1—w =constant). 
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NOTE ON FLOW IN CANALS* 
By W. H. JURNEY (Bureau of Reclamation, Denver, Colorado) 


Introduction. The design of irrigation systems to service a certain area sometimes 
necessitates the installation of pumping stations between each of which is a canal 
section with water being pumped in at one end and out at the other at constant rates. 
Starting or stopping the pumps produces long waves which are sometimes referred to 
as “surge waves” or “bore waves,” with associated changes in height of the water 
surface. Such waves have been considered by Massé! and Deymié.? It is the prediction 
of these changes in height of the water surface, under constant inflow, constant out- 
flow, or any combination of the two, which is the problem solved in this note. This is 
done by first considering an infinite canal with one source of constant inflow and using 
a method of “image” sources to produce the effect of reflections at the ends of a finite 
canal. Simplifying assumptions’ are as follows: ; 

1). The canal is of constant cross section throughout its length. 

2). The effect of fluid velocity in the canal on the wave velocity is negligible. 

3). Vertical accelerations of water are negligible. 

4). The frictional resistance to flow is proportional to the velocity. 

5). The height of water surface above normal depth, due to wave action, is small 
compared to normal depth. 

Fundamental wave equations. The usual equations of motion and continuity are:* 

On Ou B Ou On 
Ox g Ox hat 


respectively, where x is the horizontal distance along the canal (ft.), t is the time (sec.), 
h is the normal depth of the water in the canal (ft.), 7 is the height of the water surface 
above normal (ft.), « is the horizontal velocity of the fluid (ft./sec.), g is the accelera- 
tion of gravity (ft./sec.?) and B is a constant to be determined. In the computations 
on a special case for a finite canal, B was taken as igA/Q where A is the area of a 
vertical cross section of the canal (sq. ft.) and 7 is the slope of the canal computed 
from Chezy’s formula to give a flow of Q cu. ft. per sec. One alternative procedure is 
to find the frictional force resisting flow in the form* X =cu? and approximate this 


* Received April 18, 1944. This investigation was carried out by the writer in the Technical Engineer- 
ing Section of the Bureau of Reclamation which is under the direction of Senior Engineers R. E. Glover 
and Ivan’ E. Houk. All designs and investigations of the Bureau are conducted under the direction of 
J. L. Savage, Chief Designing Engineer. All engineering and construction work is under the direction of 
S. O. Harper, Chief Engineer, with headquarters in Denver, Colorado; and all activities of the Bureau are 
under the direction of H. W. Bashore, Commissioner, with headquarters in Washington, D. C. 

1P. Massé, Hydrodynamique fluviale, régimes variables, Hermann, Paris, 1935. 

2 Ph. Deymié, Proc. 5th Int. Congress Appl. Mech., Cambridge, Mass., 1938, Wiley, 1939, pp. 537-543. 

3’ Lamb, Hydrodynamics (6th ed.), Cambridge University Press, 1932, p. 254, et seq. 

* Lamb, loc. cit. 
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force by 0.75 cuQ/A. This gives the “least square fit” to the force X in the range of 
velocities 0Su<Q/A. 

It is easily seen that, for the infinite canal of uniform cross section with constant 
inflow at a point, the boundary condition is u(0, t)=Q/A =constant, where Q is the 
constant inflow. 

In order to reduce the problem to dimensionless form, it is convenient to make the 
substitutions 


B px 
p=—) pt=a, —=8, (3) 
2 a 
u(x,t) a/p) 
u(0, 0) uo = U(B, a@), (4) 
n(x, t) n(Ba/p, a/p) 
= = H 4 " 5 
(0, 0) no 
where a is the wave velocity (= V/gh ft./sec.) and 
uo = Vg/h no. (6) 
Equations (1) and (2) become 
OH ke 0U (7): OH 8) 
aB oda ( 


A solution of Eqs. (7) and (8), adapted from Heaviside’s work in electromagnetic 
theory, is 


U(B, a) = — (9) 
3 To(a 
a) = Alla) + += + (@)} 
1-3 B5{I2(a) + Is(a)} 


where I,(a) is the modified Bessel function of order n. This solution is fundamental 
to the following and corresponds to the case of an infinite canal with a barrier at the 
origin initially such that n(x, 0) =2no=const. for x $0 and n(x, 0) =0 for x>0. Upon 
removal of the barrier at ¢=0, 7 immediately drops to mo and remains there. It is this 
property of H(8, a) that makes it valuable in studying other flow conditions with 
arbitrary velocity at the origin, that is, «(0, t). In particular, to complete the solution 
for the infinite canal problem as proposed, it is merely necessary to choose 
H(0, a)=F(a) so that U(0, w)=1, and then find H(8, a) or U(8, a) correspond- 
ing to H(0, a) by an integration. 
Now an increment AF(a) put in at a=£ for 8=0 produces an increment AU at 
B=0, a=a of 
AU (0, a) & — £). (11) 


5 O. Heaviside, Electromagnetic theory, vol. 2, The Electrician Printing and Publishing Co., Ltd., 
London, 1899, p. 303 et seq. 
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Hence, after integrating and adding the initial effect, we have 


U(0, a) = f o(a — + = 1. (12) 
0 
The solution of this integral equation is easily shown to be® 
F'(a) = e€-*{Io(a) + Ii(a)}. (13) 
The effect of applying increments of H(0, a) may be represented approximately by 
AU(B, a) & (a — — (dk, (14) 


5! (a — §)? 
Thus, integration and the addition of the initial effect yields the solution of the prob- 
lem as 


(15) 


e*U (8, a) = f To{ — &)? — B} { Jo(t) + } dt + Va? — (16) 


{Tila — + — 8)} 
3! a-—é 


a—p 


1-3 Tila — I:(a — 


B* {Ii(a) 1-3 + 


3! a 5! a? 


+ e* — B{Io(a) + + - (17) 


Equation (17), on which attention is now focused, may be transformed into 
slightly better form for computation, as follows. We write 


és - 


{To(t) + 11(€) } dé 


> {Iola — s) + Ii(a — s)}ds =G,(8, a). (18) 


Therefore, 


H(B, a) = 1+ f I(t) + dt — B{Go(B, a) + Ko(a)} 


3 1-385 
{G,(B, + Ki(a)} — {G2(8, a) + Ke(a)} +--+, (19) 


6H. T. Davis, A survey of methods for the inversion of integrals of Volterra type, Indiana University 
Studies, Nos. 76 and 77, 1927, p. 51. 
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where 


{In(a) + } 


a” 


K,(a) = & (20) 


H(8, a) as given by Eq. (19) is tabulated in Table 1 for the range of values 
0 <8 SaX1, three terms of each series involved in Eq. (19) giving the results to three 
decimal places. Tables of the modified Bessel function,’ in conjunction with a Simp- 
son's rule for-five intervals, were used in making the computations. The results were 
checked by graphical integration. It is to be noted that a horizontal row in Table 1 
gives the history of the height above normal at a fixed time, while a vertical column 
gives the height history at a fixed point. 


TABLE 1. H(8, a) 


| 0.4 | 0.5 0.6 0.7 | 0.8 | 0.9 1.0 


° 
° 
° 
° 
N 
we 


000 | 
098 | 0.905 

191 | 0.998 | 0.819 

280 | 1.086 | 0.907 | 0.741 


0.991 | 0.824 | 0.670 
1.072 | 0.904 | 0.749 | 0.606 
1.331 | 1.150 | 0.981 | 0.825 | 0.681 | 0.549 
1.407 | 1.225 | 1.056 | 0.898 | 0.753 | 0.619 | 0.496 

1.480 | 1.208 | 1.127 | 0.969 | 0.822 | 0.687 | 0.562 | 0.449 

1 


.550 | 1.368 | 1.197 | 1.037 | 0.889 | 0.752 | 0.626 | 0.511 | 0.406 
.619 | 1.436 | 1.264 | 1.104 | 0.954 | 0.816 | 0.688 | 0.571 | 0.464 | 0.366 


745 
| 


a 

N 

w 


Note: The table was computed to more figures and in the range 0 $8 Sa S2, but since it is used here 
for illustrative purposes only, the abbreviated form is given. 


Finite canal section. The flow condition in a finite canal section of length Z with 
constant flow Q at one end (origin) is simulated by considering an infinite canal with 
sources of constant inflow Q located at points 0, +2Z, +4Z, ---. In the type of in- 
vestigation for which this problem was solved, the maximum height in the canal was 
the prime consideration, and only a few reflections are needed to determine this maxi- 
mum. Table 1 suffices to carry out the necessary computations. For the outflow case 
it is merely necessary to reverse results for inflow. 

Remarks. Strictly speaking, the results of course apply only to a canal initially 
at rest, and if in the case of a finite canal with a sloping bottom it is desired to com- 
pute heights subsequent to a shut down, these should be referred to the surface in 
running position which, in the case of a properly designed canal, will be parallel to 
the bottom. 

Heaviside solved the analogous electromagnetic problem for the infinite telegraph 
cable by use of operational calculus. But so far as the writer knows, the solution as 
applied to canals is unavailable elsewhere. It is seen that by the elimination of U or H 
from Eqs. (7) and (8), the differential equation which either one satisfies is 


7 Gray and Mathews, Treatise on Bessel functions, Macmillan and Co., London, 1922. 
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206 


da? | da (21) 


Thus the problem for the infinite canal, of which Eq. (19) is the solution, is equivalent 
to that of solving Eq. (21) (a, 820) with the boundary conditions 


¢(B, 0) = 0, $(0, a) = e~*[Ba{ 37o(a) + + I2(a)} + + 11(a)]. 


NOTE ON THE ELLIPTIC WING* 
By F. STEINHARDT (Columbia University) 


The acceleration potential method of Prandtl for studying the aerodynamics of 
a lifting surface is well known.! In the applications to specific surfaces one is naturally 
led to wings with circular or elliptic plan view. Prandtl’s theory has been completely 
elaborated for these two cases by W. Kinner? and K. Krienes.* 

In an effort to extend the class of wings for which numerical results have been ob- 
tained, a theory based on the work of Krienes was developed for the semi-elliptic 
wing by E. R. Lorch. Computations carried out for this case by the author had to be 
abandoned due to very poor convergence. This has suggested that the question of 
convergence in Krienes’ work, which plays an important role there, be examined more 
closely. This matter is investigated in the present note. 

We base our discussion on Krienes’ paper. The reader is referred to it for the de- 
tails which it is impossible to reproduce here. In this work the pressure p is expressed 
in terms of a potential function y. We have 


p Po poV z), Ay 0, (1) 
where p,, is the pressure at infinity, po is the density, and V is the velocity of the wing. 


In turn y is expanded in a series 


2n+1 


v=> Hav. (2) 


n=1 m=1 n=1 


(Krienes, Eq. (81) with D, =0) where 


Vn (0, ¥) 


®,, 


mm ™m m —n d n m 
bn br vn (x, ¥, 2, €)], 
m c 


* Received June 5, 1944. 

1L. Prandtl, Beitrag zur Theorie der tragenden Flache, Zeit. f. angew. Math. u. Mech. 16, 360-361 
(1936). 

2 W. Kinner, Die kreisfirmige Tragflache auf potentialtheoretischer Grundlage, Ing.-Archiv, 8, 47-80 
(1937). 

3K. Krienes, Die elliptische Tragflache auf potentialtheoretischer Grundlage, Zeit. f. angew. Math. u. 
Mech. 20, 65-88 (1940). 
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E,, being Lamé functions and 2c, 2c(1 —k?)"? beingthe lengths of the axes of the fun- 
damental ellipse. The ®, are needed to give infinite lift density on the leading edge. 
Let z=2(x, y) be the equation of the wing surface. We consider here the special 

case 
2= — x = — xao, (3) 


which represents a plane wing with a small angle of attack a». The boundary condi- 


tion is 

6z w(x, s 

bx V 


where w(x, y) is the downwash. This in conjunction with (3) leads to the following 
equations for the coefficients ay and C, in (2): 


a, = 0, (5) 
= (1 — 2 — Fa, (a=1), 
2 Cor r(2r + 1) — a(2a — 1) { 0 (a@>1), 
+ = 0, (a = 1,2,---), 


r=1 


where J,,; are certain definite integrals which are either elementary or elliptic. 

Only ®; contributes to the lift and only ®, to the pitching moment, so that we are 
interested in obtaining C; and Cy. In Krienes’ paper the series (6) are broken off 
after a=2, which leads to four linear equations giving approximations to C,, «+ +, C4; 
the results are, for an axes ratio (1—?)/?=0.2: 


Lift = 4 X 4.55aopoV? X area of ellipse, 
Pitching moment = — 1.98ac(1 — k*)"/?(po/2) V? X area of ellipse (7) 
Center of pressure at 28.3% of the maximum wing chord. 


These results seem to be in close agreement with those found experimentally. 

As a check on convergence, the computations leading to the results (7) were now 
extended by carrying the series (6) through a=3, which gives six linear equations 
in Ci, - + + , Cs. This addition of one more term leads to equations of the sixth degree 
for the coefficients of the Lamé functions (as compared with quadratics for the com- 
putations of Krienes); also, the elliptic integrals for this case give more difficulty. 
To take a>3 in (6) would make the computations almost prohibitive. 

The results now become 


Lift = 4 & 4.54a0p0V? X area of ellipse, 
Pitching moment = — 1.98ac(1 — k?)"/?(po/2) V? X area of ellipse, (8) 
Center of pressure at 28.4% of the maximum wing chord. 


This agrees closely with (7) (order of 3%). Thus one is led to believe that convergence 
of the system (6) is very rapid. 
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NOTE ON THE PROBLEM OF HEAT CONDUCTION 
IN A SEMI-INFINITE HOLLOW CYLINDER* 


By ARNOLD N. LOWAN (Math. Tables Project, Nat. Bureau of Standards) 


In a recent article! C. J. Tranter determined the heat conduction in a semi- 
infinite cylinder, in the non-steady case, by means of a combination of a Fourier 
transform and a Laplace transform. Tranter’s problem, and generalizations of it in- 
volving more complicated boundary conditions, may be solved by a method which 
was employed in an earlier paper? and involves one Laplace transform only. 

Let us consider the following generalization of Tranter’s problem: 


or? r ot 

T(r,2;4) =0 for (2) 
T(r,2;4) =0 for r=a,t>0,2> 0; (3) 
T(r, 2;t) = o(z,4) for r=b,¢>0; (4) 
\T(r,2;4) = 0 for ¢=0,aSrb;2=0. (5) 


The difference between A and Tranter’s problem lies in the fact that in A the 
boundary condition (4) involves a variable temperature. 
Let us write 


L{ T(r, 2; t)} = | e~*'T (r, 2; t)dt = T*(r, 2; p), (6) 


0 


If the system A is acted upon by the Laplace operator L defined in (6) and (7), it is 
readily seen* that the Laplace transform 7*(r, z; p) of the unknown temperature 
T(r, z; ) must satisfy the system 


(= - (8) 
or? r Or a2? k 
T*(r,2;p) = 0 for andfor r=a, (9) 
T*(r,2; p) = o*(z; p) for r= 5. (10) 
If we write 
Io(Ar)Ko(Aa) — Io(da) Ko(Ar) (11) 


— 


where J», Ky denote Bessel functions, it is easily seen that F*(a, 6,7; X) sin wz is a solu- 
tion of the differential equation (8), satisfying the boundary conditions (9) when 


A? — = p/k. (12) 


* Received June 20, 1944. 

1C. J. Tranter, Phil. Mag. (7) 35, 102-105 (1944). 

2 A. N. Lowan, Phil. Mag. (7) 24, 410-424 (1937). This article will be referred toas ANL. 
3 For details, see ANL. 
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In view of the identity 
¢*(2; p) = =f" p) sin wz sin padyda, 
and since F*(a, b, b; \) =1, for r=b it follows that the expression 
T*(r, 2; p) = = f “fF b, r; \)e*(a; p) sin ws sin waduda (13) 
0 0 


is a solution of the system A*. It now remains to subject (13) to the inverse Laplace 
operator. If 


f e~?'F,(a, b, r; t)dt = F*(a, b, r; = F* («, b,r; + *), (14) 
0 
then by Borel’s theorem‘ the inversion of (14) leads to 
2 t 
T(r, 2;2) = =f f sin wz sin f F,(a, b, r; n)e(a, t — n)dn. (15) 
rJo Jo 0 


If 
t . 
0 


the solution (15) may be written in the alternative form 

T(r, 2;2) = =f f sin wz sin padude f —G,(a, b, r; n)e(a; t — n)dn. (17) 

7/0 0 o On 

From (16) it follows that 

1 
L{G,(a, b, t)} = f e~”'G, (a, b, r; t)dt = {F,(a, b, r; t)} = b, r; 2), 
0 


whence in view of (11), 
To(Ar) Ko(Aa) Io(Xa) Ko(Ar) Y(p) 


— pt = 18 

The formal inversion of (18) is 
G,(a, br; = Aw, (19) 


Z(0) pildZ/dp) 


where the summz tion extends over the roots of Z(p) =0. Making use of (12) and re- 
membering that 


Io(x) = Jo(xi),  Ko(x) = — Yo(xi) — — + log 2)Jo(iz), 


where Jo, Yo are Bessel functions, we obtain from (19), after some simple transforma- 
tions, 


4 See ANL, p. 413. 


4 
# 
xe 


350 NOTES [Vol. II, No. 4 


em J o(a:b) 
G,(a, b, 7; = F*(a,b, r; 4) — as 
i=1 ait — Jo(aid) 
The complete solution of the system A is given by (17) in conjunction with (20). 
In the particular case where ¢(z, ¢) is a function of z only, (16) becomes 


T(r, 3; 4) = =f f G,(a, 6, r; t)e(@) sin wz sin paduda, 


which is in agreement with Tranter’s solution. 
In a similar manner it is possible to treat the more complicated case where the 


boundary condition (4) is replaced by 
(« 5:0) for r=b. 
or 


The formal solution corresponding to this boundary condition is in fact given once 
more by (17), with the function G,(a, b, r; t) satisfying the integral equation 
(p) 


(a, b, r; = (21) 
J. Z(p) 


Y(p) = Io(Ar)Ko(Aa) — Io(ha)Ko(Ar), 
Z(p) = ad{ Ko(da)Ié (Xb) — Io(da)Ké (dd) 
+ B{Io(db)Ko(ha) — 


where 


The inversion of (21) proceeds in accordance with formula (19). 


EFFECT OF A SMALL HOLE ON THE STRESSES 
IN A UNIFORMLY LOADED PLATE* 


By VLADIMIR MORKOVIN (Bell Aircraft Corporation) 


In a paper of the same title Martin Greenspan recently' determined the stress 
distribution in a large, uniformly loaded plate weakened by a small hole of an ap- 
proximately ovaloid shape. Greenspan employed a rather laborious method of piecing 
together particular solutions of the biharmonic equation for the stress function until 
all the boundary conditions could be satisfied. This process would become prohibitive 
in case of more complicated boundary conditions. It is the purpose of this note to 
apply to the same problem the elegant and more general, yet not well known, method 
for solving plane problems of elasticity which is most often associated with the name 
of N. I. Mushelisvili.? 


* Received July 5, 1944. 

1 This Quarterly, 2, 60-71 (1944). 

2 See for instance N. I. Mushelisvili, Math. Annalen, 107, 282-312 (1932). For detailed English ex- 
position see I. S. Sokolnikoff’s Mathematical theory of elasticity, (mimeographed lecture notes, Brown 
University, 1941), pp. 243-318. 
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This method rests on the representation® of the general biharmonic stress function 
U in terms of two analytic functions of the complex variable z=x-+iy, o(z) and y(z), 


U = Re f “visas (1) 


As a consequence of (1), any quantity characterizing the state of stress can be ex- 
pressed in terms of and ¥(z). Thus: 


oy + = 2[9'(z) + = 4Rel¢'(2)}, (2) 

oy — + = + (3) 

where the stresses bear the usual designation. Consider an arc PQ of some curve in 
the plane of the plate, e.g., the boundary curve, and let X,ds and Y,ds be the x- and 


y-components of the force acting on the element ds of this arc (from a direction lagging 
behind the positive tangent by some positive angle ¢ 180°). Then the quantity 


Q 


represents the resultant force on the arc PQ. When forces are prescribed on a bound- 
ary curve, equation (4) expresses the boundary condition in terms of the functions 


and ¥(z). 
In the case at hand, the equation of the boundary T is 


x= pcosB+rcos 38, y = qsin — r sin 38, (5) 


and the forces applied to it simply vanish so that the left-hand side of (4) becomes 
zero. 

It is more convenient to deal with the boundary conditions after mapping the re- 
gion exterior to the ovaloid T in the z plane conformally into the region exterior to a 
circle y of unit radius in the ¢ plane by means of the mapping function 


= w(t) = st + + (6) 


where s=(p+q)/2 and t=(p—q)/2. The values of £ on y shall be denoted by «a. It 
is noted that y corresponds to I and that ¢=1/e. Greenspan's curvilinear coordinate 
lines are obtained by setting lt] and arg ¢ equal to constants. For the sake of simplic- 
ity, the notation ¥(z) [w(¢) ] will be used. Then, 


and the boundary condition (4) reads 
+ (1/o) + = 0. (8) 


Two further results from the general theory of the method are needed before the 
functions ¢(¢) and ¥(¢) can be determined from (8). First, it is known that when there 
is no unbalanced force acting on the hole and when the stresses at infinity are uniform, 
say ¢y=Sy, Tx, the functions and take the form 


3 Derivations of formulae (1)—(4) and a formula leading to (9) are given by I. S. Sokolnikoff, loc. cit. 
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b, 
= VO + iT) +O (9) 
1 i 


where 4B=S,+.S, and 2B’=S,—S,. Secondly, one may easily derive the following 
theorems: 
(I) If f(£) is an analytic function within y, except perhaps at § =0, where it has a pole 
with a principal part of >-"Ax/t*, then 
1 As 
(II) If f(£) is an analytic function outside of y, except perhaps at § = © where it has a 
pole with a principal part of >-%A.£*, then 


1 di n 
0 


2rid, 


The actual process of solving the problem consists of multiplying equation (8) and its 
conjugate by do/27i(o—f£) and integrating around y with the aid of theorems I and 
II. The first integration yields the equation 


aq Br Bt r 
Ss 


which determines ¢(f) except for d. The coefficient a; is found by multiplying (8) 
by oda /2ri(a—f£) and repeating the integration around ‘y. One is led to the follow- 
ing equations: 


— + + a, — Br/f? = 0, (10b) 
+s B's? 4 (11) 
r—s r—s r+s 
= sBE + a,/f — Br/f*. (12) 


The integration of the conjugate of equation (8) furnishes the form of the remaining 
unknown function: 


a,=8 


— aye? + sB 


sf‘ — — 3r 


+ s(B’ + iTy)f — + Brot + + “yr. (13) 


Whenever, in plane problems of elasticity, the mapping function w({) is rational, 
the solution can be carried out in a manner similar to that above. 

It remains to verify Greenspan's results. According to equations (2), (7), and (12), 
one has 


(14) 


—av?+3rB af? + 
ort oy =2 
st4 — — 3r st* — — 3r 


which readily reduces to equation (26) of Greenspan when ¢ is on y. The equality 
of stresses at the hole and at infinity identifies the two solutions completely. 
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ON THE TREATMENT OF DISCONTINUITIES IN 
BEAM DEFLECTION PROBLEMS* 
By B.S. CAIN (General Electric Company, Erie, Pa.) 


In Mr. C. L. Brown's note on the treatment of discontinuities in beam deflection 
problems (Quarterly of Applied Mathematics, 1, 349-351), the last term of Eq. (7) is 
written in the form 


Before Eq. (4) can be applied, however, this term must be put in the form 
— PH,(x — a) — PH,(a — 3). 


This allows Eq. (7) to be written in the form 


er, 2 = — PH,(x — b) a(x — a) — a)H, 
dx? 2 2 2 
and gives, in place of Eq. (8), 
12 6 2 6 2 
Then 


+ (2b — a)§], C2 = 0. 
12d 


* Received July 12, 1944. 
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BOOK REVIEWS 


Methods of advanced calculus. By Philip Franklin. McGraw-Hill Book Company, 
Inc., New York and London, 1944. xii+486 pp. $4.50. 


By a careful selection of material the author of this handy book has put into a small volume a wealth 
of material which should enable a student to acquire some skill in the type of mathematics which he will 
want to use in his scientific work. Conformal representation, special functions and the partial differential 
equations of mathematical physics are among the topics that are treated. 

Tests are given for the convergence of integrals but the convergence of infinite series is not discussed 
at length as in many books on advanced calculus. It is thought that the student should get his knowledge 
of this subject from another text. Books have indeed been written on this one topic. 

There are useful chapters on vector analysis, complex variables, elliptic integrals and differential 
equations. Problems in maxima and minima and in the calculus of variations are also considered. 

The skill acquired by solving some of the numerous examples included in the book will give the 
would be researcher a start in mathematical analysis but he will generally find that his knowledge of the 
subject is not sufficient to enable him to complete the solution of a problem occurring in practise. He may 
need the solution of a transcendental equation or some values of a special function. The provision of 
adequate tables of functions may eventually make it unnecessary for many men to become familiar with 
the numerous properties of particular functions and the many special devices ior solving transcendental 
equations. Newton’s method is, however, explained in the text and is shown to be a special case of the 


method of the reversion of series. 
H. BATEMAN 


Modern operational mathematics in engineering. By Ruel V. Churchill. McGraw- 
Hill Book Company, Inc., New York and London, 1944. X +306 pp. $3.50. 


This is a textbook on applications of the Laplace transformation. In the first chapters the transfor- 
mation is introduced and its elementary properties are deduced and applied. In the third chapter the 
applications are made to physical problems involving ordinary differential equations, while in the fourth 
chapter they are made to such problems as the vibrating string and heat transfer (one-dimensional) in 
solids, which involve partial differential equations. 

These chapters are intended for the use of students having little or no experience with partial differ- 
ential equations. The exposition is straightforward and clear, the notation is adequate and simple, and 
the numerous problems appear to be well-chosen. Care is taken in every deduction to indicate the restric- 
tions imposed on the functions as to continuity and asymptotic behavior. An interesting section showing 
how the Laplace transformation may be used to aid in evaluating definite integrals involving parameters 
is included. Of particular importance is a section at the close of Chapter IV in which the generality of 
the methods of that chapter is inspected. The necessity for the limitation to linear equations is pointed 
out; some of the other limitations and also advantages of the methods are explained. 

The fifth chapter consists of a rather concise synopsis of the fundamental properties of functions of a 
complex variable. The various methods of recognizing an analytic function are shown; various types of 
singularities are defined and illustrated; the fundamental ideas of residue theory are introduced. This 
chapter is carefully constructed and will serve adequately for reference or review, especially for readers 
who are already acquainted with the subject. With this background, in the following chapter the contour 
integral that is the inversion of the Laplace transformation is set up and discussed. A series of ten theo- 
rems is proved, having to do with the properties of Laplace transform in the complex plane and the valid- 
ity of the inversion integral. The restrictions on the function and transform in most of these theorems 
are rather severe; it seems possible that this is due to the author's having attempted to achieve general- 
ity without using advanced mathematical methods. Useful methods of deforming the contour of inte- 
gration in special cases are illustrated. 

In the next two chapters the complete theory now available is applied to boundary-value problems 
in heat conduction and mechanical vibrations. These applications are characterized by rather extreme 
attention to the restrictions imposed in the theorems of the preceding chapter. The author points out 
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that the solutions obtained by formal use of the transformation and its inversion can be verified by 
showing that they satisfy the conditions of the boundary-value problem; nevertheless, in some examples 
he analyzes them at great length to show that they satisfy the conditions required at infinity, etc., for the 
theorems, and therefore are valid solutions. It seems doubtful that this material is useful in an engineering 
course, since engineers are unlikely to adopt such techniques in preference to the more direct procedures. 

The reviewer notes the absence of any explicit statement of the well-known methods of determining 
the initial and asymptotic behaviors of a function by expanding its transform for large and small values 
of its argument, respectively (see, for example, Doetsch: Theorie und Anwendung der Laplace Trans- 
formation; Berlin (1937) p. 243). In engineering applications these methods have been found to be ex- 
tremely useful. 

The final chapters are concerned with Sturm-Liouville systems and finite Fourier transforms, respec- 
tively, although these subjects are not intimately connected with the Laplace transformation. This 
discussion of Sturm-Liouville systems, while concise, represents a sort of generalization that engineers 
and physicists will find instructive and useful. The finite Fourier transform, however, seems to the re- 
viewer to be a rather artificial conception of little practical importance. 

The textbook includes tables of finite sine and cosine transforms, a table of operations (properties 


of the Laplace transformation), and a rather complete table of Laplace transforms (122 entries). 
W. R. SEARS 


Chemical engineering nomographs. By Dale S. Davis. McGraw-Hill Book Company, 
Inc., New York and London, 1944. ix+311 pp. $3.50. 


This book is a collection of over 200 nomographs relating to problems in chemical engineering. The 
accompanying text consists of a discussion of the validity, application, and limitations of these charts, 
and directions for their use. As a reference textbook for chemical engineering coursework, it would be 
particularly valuable as a companion volume to the earlier work, “Empirical Equations and Nomog- 
raphy,” by the same author. To a practicing engineer, the most important feature of the book is the 
excellent selection of nomographs. This selection was made from a total of 600 sources; “in most in- 
stances, the hundreds of requests for reprints . . . during the past 14 years have guided the choice.” 

The 19 chapters of the book may be divided roughly into three types: 

(1) Those dealing with unit operations (e.g., flow of fluids, distillation, extraction, etc.). 

(2) Those giving a specific property (e.g., density, viscosity, thermal conductivity or vapor pressure) 
of a wide range of substances. Vapor pressure-temperature-concentration relations are particularly well 
covered. 

(3) Those giving a wide range of properties for a particular set of materials; e.g., acid nomographs, 
milk and cream nomographs, and a large group devoted to problems of the paper industry. 


The charts are well grouped and indexed and are printed in a form large enough to be read easily. 
R. LEININGER 


Tables of Lagrangian interpolation coefficients. Prepared by the Mathematical Tables 
Project, Work Projects Administration of the Federal Works Agency; conducted 
under the sponsorship of the National Bureau of Standards. Official Sponsor: 
Lyman J. Briggs. Technical Director: Arnold N. Lowan. Columbia University 
Press, New York, 1944. XXXVI+392 pp. $5.00. 


The main tables give the Lagrangian interpolation coefficients for interpolation polynomials of the 
degrees two to ten. The following list gives the range and interval of the argument, and the number of 
decimal places in the usual “Sone: the arguments of the equidistant given values of the function 
being denoted by 0, +1, +2,--+-: three point formula [—1(.0001)1; exact], four point formula 
[ —1(.001)0(.0001)1(.001)2; 10D], five point formula [—2(.001)2; 10D], six point formula 
[- —2(.01)0(.001)1(.01)3; 10D], seven point formula [—3(.01) —1(.001)1(.01)3; 10D], eight point formula 
[ —3(.1)0(.001)1(.1)4; 10D], nine point formula [ —4(.1)4; 10D], ten point formula [ —4(.1)5; 10D], eleven 
point formula [—5(.1)5; 10D]. Additional tables give three- to eight-point interpolation coefficients for 
intervals of .1 particularly suited for the purposes of subtabulation, three- to eight-point interpolation 
coefficients in fractional form for intervals of 1/12, and Lagrangian integration polynomials and co- 


efficients. 


W. PRAGER 
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A treatise on the theory of Bessel functions. By G. N. Watson. Second Edition. 
Cambridge University Press, Cambridge, England. The Macmillan Company, 
New York, 1944. viii+804. $15.00. 


For those already familiar with Watson’s “Bessel Functions” a new edition needs no recommenda- 
tion, and indeed the fact that a second edition has appeared in the midst of wartime difficulties is in itself 
sufficient evidence of the book’s importance. Only a few minor corrections have been made at this time, 
but it is still cause for rejoicing that “Watson” is once more readily available. In the preface to the 
new edition the author regrets that his interest in Bessel functions has waned since the book was first 
published, and rightly remarks that to have included all the new material which has appeared in the last 
twenty years would have meant rewriting most of the book. It is unfortunate that he did not find it 
possible to bring at least the bibliography up-to-date so that for recent work we must still look elsewhere. 

Nonetheless, the book remains an invaluable compendium of information for any mathematician 
whose work ever touches on Bessel functions. While the presentation is mostly designed for the pure 
mathematician the applications are not completely neglected, and the applied mathematician will find it 
particularly useful if he regards it as a treatise on the theory of functions of a complex variable, with 
applications to Bessel functions. 

If there are any applied mathematicians whose work has not so far brought them into contact with 
“Watson” a brief list of some of the topics treated may prove illuminating. Such a list, by no means com- 
plete, would include: solutions of Riccati’s equation, expansion of functions in series of Bessel functions 
(including Fourier-Bessel, Dini, Kapteyn and Schlémilch series), addition theorems, methods of evalu- 
ating definite and infinite integrals, asymptotic expansions for Bessel functions of large argument and of 
large order (using the principle of stationary phase and the method of steepest descents), determination 
of the zeros of Bessel functions, and last but by no means least, tables of various Bessel and related 


functions. 
Both author and publisher are to be congratulated on the successful reissue of this classic treatise. 


May it long remain in print. 


M. C. Gray 
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By T. Wittrraker. 4th edition, 460 
pages, 6 x 9, Originally published at 


$6.00 $3.50 
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